System and method for segmentation and automated analysis of the three-dimensional structure of duct

专利摘要:
System and method to quantify, based on a computed tomography, the metric and morphological characteristics of a conduit of irregular shape and heterogeneous walls and content. It is characterized by automating the quantification with the following stages. Dimension in each image a duct environment; segment it automatically iterating the triple process of extraction of traits (414), classify the pixels (416) and rectify the classification, adjusting it to the existing knowledge about the properties of tomographed objects (418); the iteration stops when the segmentation converges (420); the segmentation procedure is repeated with the other sections of the tomography (424) to obtain volume segmentation (426). With it and the tomographic metadata, the models of the segmented conduit and its axis are automatically created in the physical coordinate system; and from these the measurements and morphometric descriptors that characterize the three-dimensional structure of the conduit are automatically generated. (Machine-translation by Google Translate, not legally binding)
公开号:ES2608037A1
申请号:ES201500707
申请日:2015-10-01
公开日:2017-04-05
发明作者:Lucia Jañez Garcia
申请人:Lucia JÁÑEZ GARCÍA；
IPC主号:

专利说明:

TECHNICAL SECTOR
Quantification of the dimensions and shape: sectional and three-dimensional ducts through automated amilysis of your CT images.
BACKGROUND OF THE INVENTION
The technical problem is focused on the dimensional and morphometric quantification of irregularly shaped three-dimensional ducts and of heterogeneous walls and content, as well as their sections, by means of the automated analysis of their computed tomography images affected by noise and generated by any of the following modalities: computed tomography (CT) and microtomography (Il-CT); conical beam computerized lomography (CBCT); magnetic resonance imaging (M R 1) structural, functional and diffusion; positron emission tomography (PET) and combined studies, of this with other modalities (PET-CT, PET-MRI); magnetoencephalography (MEG); optical coherence tomography (OCT) in the temporal, spatial and Fourier domains; focal and fluorescence microscopy with one and two photons; tomosynthesis and other lomographic modalities in use and in development.
In living beings there are several types of irregularly shaped ducts and of non-homogeneous content and walls to which this invitation applies in its different embodiments and that include the nasolacrimal ducts, mandibular duct, paranasal pits, ducts
auditory and other ducts.
The term 'oconduct' is used here to restrictively refer to conduits characterized by their irregular fomla and for having heterogeneous walls and contents. But regardless of the nature of the material that delimits it (bone, living tissues, rock) and the content inside that can be solid, liquid, g, aseptic or mixed and can be totally or partially in motion or static (such as air, tear fluid, solid and mucous deposits in a tear canal or the fossil of a duct embedded in the inside of a heterogeneous rock).
There are already general methods that have proven their usefulness to segment the sections of a conduit of approximately constant section and homogeneous walls in a tomography. such as blood vessels in a computerized axial tomography for clinical use: the thresholding of gray levels, the detection of contours based on the first and second partial derivatives (Jang et al. 2014), the active contours (also called amadossnakes) (Kroon. 2011). the growth of regions (Kang el al 2003), the methods of mathematical morphology (erosion, dilatation, openness and closure. with structural elements of different sources). the "basin" ("watershed") methods. those based on a ", texture analysis or level sets. the methods of classification, are some of the most used; but when we have applied them to segment small and complex ducts, specifically human lacrimal naso ducts in axial tomography images computerized for clinical use, each method in itself has proved incapable of making a segm <: correct and robust automatic nestation; in fact we know it has not yet been published until now any procedure to automatically segment them. neither can they make any of the existing methods viable on their own; specifically those based on Gaussian filtering solve some problems (such as variability due to noise) but complicate others (extending the very low contrast sections in the contours of the duct): Anisotropic filtering has a disturbing effect due to the importance of noise in tomographic images Thus, highlighting contour lines due to noise and eliminating pieces of real contours due to their low contrast.
The difficulty with ducts such as tear ducts comes from the combination of the following factors: the small size that the tear canal section can have in relation to the size of the image, which is sometimes around 10% of the number of rows or columns from image; the variability of the gray levels within the channel: the channel itself - especially when it is clogged and in the "walls that delimit: the great fluctuation of the ConTera along the sectional sectional channel that in some areas can even reach be canceled: the low spatial resolution of the tomography - around 0.3 mm - in relation to the thinness of the bone wall that in some areas delimits the canal and that only in some areas only reaches a thickness dc O '1 mm; the random fluctuation that the values of the voxels undergo in all the tomogrMicas modalities that is linked to the signal-to-noise ratio and that generates false contours whose contrast is sometimes significantly greater than the txx in some sections of the actual contours; the shape irregular and changing section along the canal and between them.
In addition to the tear ducts cxistt ~ n in living organisms other ducts, also irregular and non-homogeneous walls whose study poses similar problems and where the situation of the technique is similar, so they also benefit from the technical progress provided by this invention and that include the mandibular duct. paranasal pits and others.
The difficulty that for the reasons described presents the segmentation of ducts with one or several of the characteristics described in the previous paragraph has forced to resort systematically to manual methods to carry out not only the segmentation but also the dimensional measurements and their shape.
Manual segmentation is still practiced today (Estes et al. 2015), and to compensate for the subjectivity introduced by each expert operator, the same boundaries have been repeated with several experts. whose number is usually in volume two or three due to the work and time: time it requires; but repeating a delimitation or measurement by experts is not the solution, because although it can compensate for random errors it does not correct the systematic errors introduced by each of the experts; In addition, repeating a delimitation or measurement such a small number of times reduces very little the confidential interval for the average value obtained, and that lack of reliability inevitably affects all measurements that directly or indirectly rely on the aforementioned manual operation and degrades the reliability and repeatability of the results.
Another problem. associated with manual segmentation, is that the dimensional measurements of them are also performed manually. measuring lengths and angles between points of imprecise location visually fixed both in the tear ducts (Yong et al 2014, Takahashi et al 2014) and in the mandibular canal (Pyun et al., 2013). measuring the distance between two points of the channel by the number of cuts that separate them (Takahashi et a12014) or by the length of the straight segment that joins them although the channel is curved (Ramey et al. 2013). using ellipses to measure the area of irregular figures that are actually closer to the ovoids. and using ellipsoids to measure volumes of cavities whose shape is irregular. Manual measurement means that the measurements require a lot of time and are carried out in fewer points than necessary and in subjectively fixed positions; In addition, the accuracy of the measurements obtained is reduced due to the motor limitations of the person who takes them manually, the subjectivity when choosing the points used in each measurement, and the discrepancy between the geometric models used and the irregular shape of the real structures; the measurements are also different [due to the inaccuracy in the visual representation of the ducts to be measured when they are too small or when their contours are blurred when enlarged or perceived displaced by the effect of inadequate ambient lighting. The methods for evaluating volumes that are based on manual segmentation of the duct in the cuts of rAe (Estes et al 2015) manage to gain precision with respect to those based on approximations by means of geometric bodies, but continue to consume so much time that hinder its use in clinical practice and also continue to suffer from the lack of precision and objectivity in the outline of the contours as a result of human motor and perceptual limitations; In addition, the problem is exacerbated when the contours that appear in the image are not well defined or too small in relation to the resolution of the tomography.
A third consequence of the difficulties cited for automatic segmentation is that the study of the form is also supported by manual operations on the image, being therefore affected by the subjectivity and lack of reliability that this implies both for the descriptors of classical form and for derivatives of geometric morphometry. In fact, regarding the description of the fan of the sections of the duct (20), the duct itself (3D) and its dynamic evolution ot: volutive over time (40), there are two approaches: the classic and the most recent based on geometric morphometry. In the classical approach to flat shapes, manual measurements of some of its dimensions are made: maximum and minimum diameters, areas, etc. (Takahashi on 2014) and intuitive indices such as eccentricity or elongation and abstract mus indices such as standard, central, normalized and invariant central moments, or Fourier descriptors and others defined in Gonzalez and Woods (2002), Levine are calculated (1985) and Hu (1962). The shape of the three-dimensional ducts is sometimes described by approximations in terms of solids on a regular basis, such as cylinders or ellipsoids, given the values of their dimensions; however, for very irregularly shaped ducts, which is common in biological structures, these classic nuances are too inaccurate and more precise descriptions of their fanna and three-dimensional structure are necessary. The geometric morphometry generates for the flat and three-dimensional forms indexes or descriptors qU (: in some cases they are based on reference points
or marks (Iandmarks) and in others in the contours sccionalescs or in the surfaces that delimit the conduit; for the first ones the analysis of marks has been developed (Bookstein, 1991) and for the contours we have used methods of adjustment of curves (Rogers, 1982), analysis of main components (Glassburn, 1995), analysis of classic Fouricr (Temple , 1992) and elliptical (Ferson et al., 1986), analysis of eigen · fonnas (Lohman, 1983), etc. But the great precision that a mathematical analysis of the form can offer is degraded when the location of the reference points has to be done manually. claiming again the automation for the analysis of the form in 20 and 3D, as well as for the automated obtaining of the data or marks that serve as a starting point for the geometric morphometry.
A fourth consequence of not having automatic segmentation methods is that of doctors. veterinarians and other clinical professionals lack tools to visualize the tear ducts; Researchers have to work with mental models of the ducts that can hardly be validated and contrasted. based on regular geometric figures such as ellipses and ellipsoids to model objects such as the nasolacrimal canal (CNL) that is very irregular and rarely ellipsoidal. and they are forced to work with three-dimensional digital models so primitive and distant from those that allow the current graphic techniques that make them scarcely interpretable and useful, as evidenced by the figure recently published by Estes el al. (2015).
The problem posed by manual segmentation in tomography does not arise only in medicine but also in other areas where tomographic images represent a tool of frequent use and increasing importance: veterinary, biology, anthropology, etc. From an economic and industrial perspective. User intervention reduces the speed of the analysis and increases the associated personnel costs. As the use of measurements on digital images is increasing, there is also a growing need to develop methods that allow more measurements to be made, in less time. more accurately and with more objectivity. Thus, finding an efficient <.] Oration that automates the segmentation process. modeling, measurement and morphometry in tomographic images. <; e has become a manifest necessity and satisfying it is an objective of great scientific and industrial interest.
This patent contributes to this by providing a method to determine the three-dimensional structure of ducts from their tomography. which allows the duct section to be automatically and objectively segmented into tomographic images, to generate quantitative models of the duct from that segmentation. and on them, carry out the dimensional and morphometric measurements that characterize the quantitative level and the objective of the tomographed duct.
EXPLANATION OF THE INVENTION
In this document the use of the Leninins "a ~," a "," some "and" some "and similar references in the context of the description of an embodiment of the invention, and especially in the context of the claims. to take to cover both the singular and the plural, unless otherwise indicated in this document.oo it is clearly contradicted by the context. The terms "includes", "comprises", "that have", "including", " understanding "and" with "should be interpreted as open terms (ie, which means" including, but not limited to, "), that is, they should not be construed as exductors of [the possibility that what is described and defined includes more elements, stages, etc. unless otherwise indicated, the mention of ranges and ranges of values in this document constitutes an abbreviated method of individually referring to each value within the range or interval unless otherwise indicated. what would be included in the eSle document, and say This value is incorporated into the specification as if it were cited individually in this document. All methods described in this document may be carried out in any suitable order unless otherwise indicated in this document or clearly contradicted by the context. The use of examples, or t: related expressions (for example, "as is"), seeks only to clarify the invention better and does not plan a limitation on its scope. Preferred embodiments of this invention represent a known way of carrying out an embodiment of the invention, but there may be modifications of these preferred embodiments that are apparent to those of ordinary skill in the art and hence the invention can be brought to the pni (; ethics differently than specifically described herein. Accordingly, this invention includes all modifications and equivalents of the subject matter in the claims appended thereto as permitted by applicable law. On the other hand. any combination of The elements described above in all possible variations of the same are encompassed by the invention unless otherwise indicated herein or clearly contradicted by the context.
This invention proposes a system, a device (fig. 1) and an automated method and method (fig. 2) to obtain in an automated way the dimensional quantification and metric morphology of irregularly shaped three-dimensional ducts and of heterogeneous walls and content, as well. as of its sections, by means of the automated analysis of its tomography images affected by noise. The method comprises the procedure of segmenting a duct in the tomography, from that segmentation create numerical models of it, and on these models make measurements of the dimensions and shape of the duct. In order to simplify the exposure
in this document we will call the physical space delimited by its walls "conduit"; the rest of the space captured in the tomography will be called "background", although it may include different types of substances or tissues.
THE METHOD
The described method (fig. 2) comprises the e "lapas described below: obtaining the lomography (200). Preprocessing (202), segmentation (204), modeling (206). And measuring (208) morphometry (110 ) And results (212).
This amilysis sequence described below is intended to analyze a single conduit, but this does not represent a loss of generality. since when the tomography contains several ducts (as is the case of the two tear ducts in a head scan) then the tomography will be processed as many times as the ducts contain. analyzing only one each time. In the case of ducts whose structure is complex. With the aim of differentiating several integral parts of the same conduit, the processing will be done considering each part as a different conduit and subsequently the resulting conduction will be taken as the junction of the segmented conduits separately, and as characteristics of the complete conduit, those resulting from the integration of the data found for its components.
1. OBTAINING THE TOMOGRAPHY The first stage (200) of the method comprises obtaining the tomography. and maximize the number of cuts when required:
1.1. Acquisition of the tomographic volumeJlL (200). The tomographic volume can be obtained directly from the equipment that performs the tomography (108) following the specific procedure of the modality used (CT, NMR, OCT, etc.); it can also be taken from an rACS or other storage system where it has been previously archived (1 06). Figure 8 shows as an example an integrale image of a human head tomography; In its upper part the arrows indicate the location of the sections of the two tear ducts. This image is digital and consists of a mosaic of 5 12 rows and 5 12 columns. where each tile of the mosaic has a homogeneous gray level taken from a scale that differentiates about 4000 different gray levels (dynamic range) :. when the image is presented in size like the one in the figure, the mosaic elements are t ~ 1O small that do not differ from each other and originates
the sensation of a continuous image: the tomographic volume is formed by a sequence
close to one hundred images stacked neatly and each of which represents the radiography of a virtual slice of the head at different equally spaced heights.
Each digital image is represented by a numerical matrix of variable dimension and dynamic range according to the tomographic mode used, and the volume of tomography (hereinafter we will simply call "tomography") is shaped by an ordered sequence of these matrices that together form a hypernatrix that can have any number of dimensions. but typically it has three spatial dimensions and sometimes a temporal one. Each element or voxel of said hypemlatrix represents the average of a property that the tomographed object possesses in the physical volume that corresponds to that element and that is usually a cube or a square-based cousin that fan part of that object.
The geometric characteristics of the tomographic acquisition that pennilcn the change of coordinates from the image space to the physical space (and vice versa) are given in the tomography mctadata: these metadata also contain ample complementary information.
1.2. Maximize the number of cuts that include the duct. Sometimes the analysis can be carried out based on the reconstruction generated directly by the tomography equipment: this is the case of the tear ducts in clinical tomography, where generally the angle between the cutting plane and the longitudinal axis of the duct moves away excessively from perpendicularity. When this is not achieved during the acquisition of the images. then a multiplanar reconstruction will be generated as an initial step that maximizes the number of eons on the duct to capture as much tomography as possible on the duct.
2. TOMOGRAPHIC VOLUME PREPROCESS (202)
(a) This stage of the method includes creating the files with the image data and metadata; normalize the figure-bottom gradient, obtain the binary mask of the region of interest (RDI) that encompasses the section of the duct in each image, and determine the values of the general parameters associated with the nature of the object and its tomography (fig. 3).
2.1. Create the files with it !; Image data and metadata (302). A tomogrMico study generates two types of data: images and metadata. The images contain the representation obtained from the part of the tomographed body: the meladatos, also called "headers" of the images. they contain information that facilitates the analysis and interpretation of the image data: number of rows and columns of each image, number of bits / pixel, physical dimensions of the voxels, width and spacing of cuts, origin of physical coordinates, sense of acquisition of the tomography. type of core used for the reconstruction of the volume, and sludge with other data specified by the data dictionary of the format used. Each image has its own header. and the tomography equipment is supplied both together in what is called a study or volume. organized according to the specifications of the fonnate used and which for medical imaging is often DlCOM. In this stage, images and metadata are separated into different files, maintaining the biunivocal correspondence between each image and its header: the set of all images is grouped in a V-Hipella with 3. 4 or more dimensions; its elements vll ~ will have the same number of indices: the first two indicate the row and the column they occupy in the image of a cut, the third represents the number of the spatial plane in which the cut was made. the fourth index (when it exists) usually indicates the order number of the moment1: Or in which the tomography was acquired. The set of
the metadata is grouped in an M matrix.
2.2. Standardization of the gradif: nte figure-bottom (04). The method of this invention assumes that the average of the pixels per1 (! Necent to the conduit is less than that of those belonging to the background in an environment close to the conduit. ESI0 does not represent a loss of generality in its application, since in the case that the tomography obtained does not comply with that assumption then it will be negativized and the metadata that inform the transformation to be applied to obtain the values of the image will be modified accordingly, which in DlCOM are in the labels (0028. 1052) AND (0028.1 053) .
2.3. Define the region of interest (RDl) in each court (306). The Region of Interest or RDI is understood here as the area of the image whose shape may be irregular and which will include the entire section of the duct and an additional part of the surrounding substance; said region will be defined by a mask that has the value 1 at the position of the pixels that belong to the region of interest and Oen the rest; the delimital;: ion can be done by automatic, manual means
or assisted by a computer. The delimitation of the region of interest constitutes a totally different problem from the segmentation of the duct and it is proposed to do it separately: in the case of images with simple structure, its realization can be trivial. defining the ROl equal to the whole image; in other cases. as in the head cuts that do not present too much structural complexity locating the ROl is a problem that can be addressed automatically; and in very complex images the automated solution can still take time and in the meantime you also have to resort to manual or assisted delimitation. In any case, the limitation of the RDI must be done in a way that does not interfere with the automatic segmentation of the duct, for which it will include a strip of the surrounding substance, which should be wide enough and at the same time as homogeneous as possible. and in such cases it is necessary to resort to a manual or assisted delimitation of the RDI. Fig. 9b shows an example of an assisted embodiment: the ROl is given by the polygonal of thick stroke; The fine lines are help lines whose purpose is to prevent the delimitation from entering areas where it could interfere with the segmentation of the conduit that will be carried out later in an automated manner.
The assisted delimitation of the regions of interest for each tomographic image includes one or more of the following aids: zoom with interpolation, smoothing of the image. superimpose reference lines on the image. and visualize the image with a level in a range of background values close to [the gray levels inside: 1 duct and with a window in volume one third of the range of background values.
The automatic delimitation includes the mart: aje of an initial point or semi-pan to panir of which by approximate techniques of growth of regions an approximate delimitation of a region that encompasses the conduit is determined; this delimitation is spreading to successive dogs based on the contours identified in these; such delimitation does not require the precision of segmentation for metric purposes; Finally, a widening of the region found in each cut is practiced to avoid interfering with the segmentation process and undergoes a verification process.
The result of delimiting the ROl in an image will be a binary matrix of the same dimension as the image and with the value 1 in bp: eles included in the RDI and Oen the rest; The resulting set of all the RDI masks obtained for all the images of the lomography will constitute a new hypermatrix R that is related to the hypermatrix of images V from which they come from the bijective correspondence g: VR ta l quef {v '/ klJ == r'JkI; for the images in which the conduit is not represented, the corresponding masks of the region of interest will give sludge their equal elements to zero.
Separating the methods to delimit the regions of interest and the segmentation of the conduit within them significantly reduces the time of calculation of the segmentation and simplifies the segmentation process.
2A. Determine the values of the 2eneraJes parameters (308). The automatic segmentation method requires that information be provided on the nature of the material to which the tomography corresponds and on some characteristics of it. This information contribution is specified in setting the values of a number of parameters that will be used during the automated segmentation process. In general, the general parameters include: a) the numbers of the first and the last image of the tomographic series that intersect the object: b) the nature of the index of similarity between two successive segmentations,
c) the value of cone required at said index to stop the segmentation iterations with each image; d) criteria to determine at runtime the form and size that the structural elements must have in the operations of mathematical morphology and that depend on the pre-existing knowledge about the shape and curvature that the surface separating the conduit can take. of its surroundings; and e) the ranges of values envisaged in the tomography for each type of substance that integrates the tomngratiado ohjeto, derived from pre-existing knowledge.
Other auxiliary data will also be entered in the “! Rales” parameter file. such as the name and address of the physical files that contain the images and their meladatos. the initial and final numbers of the eons to be processed. the names of the files where the program must record each of its results; and other auxiliary data that obviously need not be detailed. 3. SEGMENTATION (204)
To obtain the segmentation of the duct in each cone and in the volume. the method described below will be applied and that includes iterating the following sequence of methods (fig. 4):
a) Evaluate and record the specific parameters for the current iteration (4 l 2). which generally include: a) determine the reference region. that for the first iteration of segmentation in a lomographic image is ROl, and for the second and subsequent ones it is determined by the rectified duct mask obtained in the preceding iteration; b) determine the area of the reference region in the physical coordinate system, using information on the size of the ~; pixels contained in the metadata: e) calculate the physical coordinates (e, c¡) of the central reference point, given by the centroidc of the reference region: and d) determine the shape and size of the structuring elements to be applied in mathematical morphology operations - erosion and dilution - calculated according to the size of the reference region and the criteria established in the general parameters.
b) Calculate the feature vector associated with each pixel of the region of interest (ROl), a vector that will have at least dimension 2. each of whose coordinates is the resolute value of applying the corresponding feature extraction operator on said pixel, calculated from the value of the pixels and the general parameters and the specific parameters of each iteration that are recalculated at the beginning of it: (414). This step consists in calculating the feature vector associated with each pixel of the ROl by applying in each pixel the set of feature extraction operators. each of which will generate. the corresponding coordinate in the vector. using for this the general and specific parameters. The result of this stage is the matrix of features F, which is a matrix of the same dimensions as the image and whose clement / y is the vector of features corresponding to the pixel located in the same position ij of the image.
Each operator can incorporate as a preliminary method to imitate a specific noise reduction network method: this noise elimination can also be incorporated as a previous and common phase to all of them by means of a medium filter, or a Gaussian low pass filter or the convolution of the image with an optimized matrix to limit the specific type of tomography noise.
The feature vector includes as one of its coordinates the value calculated in each pixel by a non-linear, 2-D or 3D operator, whose parameters are dynamically recalculated for each segmentation iteration; and as another of its coo: the value adopted in the said pixel by a discrete differential operator calculated using the physical coordinate system determined by the metadata.
The non-linear operator used for the (traction of features in a first embodiment of the invention for tear ducts has been a threshold with two thresholds. I1 and 11. whose values simultaneously simultaneamenle the general parameters set a priori and a statistic of central tendency of the values of the image in a crazy environment l of the central reference point; said operator calculates for each pixel the value of the feature based on the value of the said pixel, xy, assigning the value 'T' when x, ) <I}, the value "O" when Xy> I} Y the value "0.5" otherwise. Thresholds 11 and tl are calculated as follows: I1 = min [Cmax. max (C, .. ,,,, MJ) J y / 1 = (Bmm T (¡) / 2. where Md is a central tendency statistic of the values of a submatrix nxn of the image centered on the central point of reference, taking n a value around 5; Brnm represents the minimum anticipated a priori value for pixels not belonging to the conduit. and Cn "" and C ", ax respectively represent the lower and upper COLas respectively for the range of the ma value .: l (imo provided in each eorte for pen pixels: in the "duct": the values of C ", ax, Cm¡ and B" ,,,, are established by an expert depending on the nature of the tomographed duct and of the tomography modality used. The parameters Cm "" Cn, Q; t. and Bmm take values in the intervals whose approximate limits are respectively 30-80, 80-120 and 300-700.
In a second realization the linear lilac operator consists of an active contour that a) is initialized with the same shape and center as the reference region and with a size equal to a fraction of it: b) the contour evolves towards the periphery of the channel attracted by a function of the partial derivatives of the pixel value until it reaches its steady state; and e) the nonlinear operator assigns the value 'or O' to the trait according to whether the pixel is located inside or outside of the stationary contour_ and 0.5 otherwise.
In another embodiment, the nonlinear operator consists of an operator based on watershed C level sets that a) applies a segmentation algorithm to the RO using the watershed level sets technique; and b) the operator assigns the value 1, 060'5 to the trait according to whether the pixel is located in the greater of the interior regions, in the greater of the peripherals, or in another. respectively.
In the first embodiment the differential feature space operator consists of a Laplacian operator that cn each pixel of the RDI assigns to the corresponding feature the discrete Laplacian value of the image values over fd cited pixcl.
In another embodiment the differential feature space operator consists of the directional derivative operator calculated by the discrete derivative in the direction of an identifiable reference point located inside the conduit.
There may be more than one non-linear operator and more than one differential.
e) Classify each pixel of the ROl by labeling it with the value l or O, representative respectively of "conduit" and "background", resulting from applying a binary classifier to the feature vector of said pixel: and creating the matrix of the initial mask of the conduit with the values thus obtained in the pixels of the ROl and with O in the others: 416). The input data for this stage are those of the array of features F. and the output is a binary matrix of the same size of the image that has I in the position of the pixels classified as "object" and O in those classified as " background ", matrix that we will call the initial duct mask.
In the first embodiment, the classifier e, s a Detentionist decision tree that labels each pixel of the ROl as a "conduit" or as a "background" based on the values of its feature vector and which comprises the following steps i) label as " duct "the pixels in which cI 5 feature vector takes the value I in the coordinate generated by the nonlinear operator; label as "background" those in which said coordinate is worth O; and label those in which it is worth 0.5 as a "conduit" or as a "fund" according to whether the coordinate from the differential operator is positive or negative, respectively; ji) classify the remainder of pixels with 0.5 in its first coordinate as a "conduit" or - background ", when in a small environment there is a majority of
10 pixels already classified as "conduit" or "background", respectively; and when both proportions are equal then it will be classified as "background" when the average value of the pixels in the aforementioned environment is closer to the higher threshold and as "conduit" in another case d1 Rectify the classification given by the initial mask of the conduit reversing the binary value of its elements corresponding to the ROl when they do not fit the knowledge
15 pre-existing on the nature of the substances constituting the tomographed material and on the shape of the contact surface between the duct and the bottom; the matrix thus obtained is the rectified duct mask; (418). At this stage, pixels of the duct and the bottom are re-classified which do not meet the conditions imposed by the physical properties of the tomographed matter; in the case of a conduit, the following conditions are imposed: a) that they do not exist
20 isolated bottom regions inside the duct, nor vice versa; b) that the curvature of the contact surface between the duct and the bottom (taken in any direction and in particular in that marked by the lomographic cuts), does not exceed the known limits and consequently there are no fibrillar penetrations of the duct on its walls or vice versa; e) that both the region classified as "conduit" and the region classified as background "are both regions
25 related and that its union includes the entire image; and d) that the conduit does not have points of contact with the set of complementary pixels of the region of interest.
On these premises the method of rectification, ation includes the following steps: 1) apply a morphological closure, consisting of reiterating, with a limit in the number of iterations, a dilation 30 followed by erosion, with surrounding structuring elements whose radius It is in inverse relation to the maximum curvature admitted in the separation surface of the · -figure "and the" bottom "and in direct relation to the size of the reference region: 2) apply a morphological closure with a limit on the number of iterations and with fibrillar and small-sized structuring elements; J) when the set of pixels classified as "conduit" forms several unconnected regions 35 then all pixels belonging to the most regions are reclassified as "background"
small, retaining the classification as "conduit" only for those pixels that belong to the region that contains the largest number of pixels classified as "conduclo", and if there were more than one region with maximum clamping, the one that has its centroid closest to that of
the region of interest; 4) reclassify as "" background "the pixels that meet the following conditions manually: i) the coordinate of their feature vector generated by the non-linear operator with the value 0.5, ii) is on the periphery of the conduit This stage of rectification by prior knowledge is repeated twice, placing the four steps in each order in a different order to minimize the interactions that result from the order in which the different ones are applied. " algorithms Control the number of iterations (420). In the iteration control stage, the current rectified duct mask is compared with the rectified masks obtained in all previous iterations therein. image; if the pre-established similarity criterion is not satisfied in any comparison. en1: Onces continue in the initial stage to do an additional iteration; otherwise add the current rectified mask. as a definitive duct mask, to those obtained previously, to build with them all the binary hyper ~ matrix that defines the segmentation of the duct. Thus, the binary duct mask is incorporated into the binary hypernatrix B that collects the duct masks in all cuts; additionally, the user is shown pma tracking and recorded in permanent and computer-readable media for later use by modeling, measurement and morphomctric methods. In the "hypem, binary atrium" 8 each binary image occupies the same position it occupied in the hypermatrix A the image of gray levels from which it comes. Consequently B will have the same structure as A, and the indices of each clement of the binary hypermatrix have the same meaning as in the matrix of the image "the first indicates the row. The second the column, the third the depth, the fourth the time, and so on; the difference is that each element of this hypermatrix can only turn two different values, which may be any, but that are generally O and 1 (the exposure follows simply under this assumption): it takes the value I when the corresponding pixel has been classified as "conduit" and takes the value O when it has been classified as 'bottom'.
During the automatic segmentation process the program will record in a log file
(130) And it will make the user accessible on a screen or by other means (128, 102 or 104) the relevant data on the segmentation obtained in each iteration for all cuts, without requiring any response from him, but allowing him to do so through the mouse or other available means of entry (126) to take control of the decisions, either by accepting the
segmentation resulting from a certain iteration or stopping the execution of the program, The aforementioned relevant data will include at least the identification of the cut, the number of the iteration, and the values of the similarity indices between the segmentation resulting from that iteration and those obtained in previous iterations,
Advantages. The proposed segmentation method is robust, it offers very concordant results with the manual method. and it surpasses it in terms of objectivity, reliability and performance: it is also deterministic, because it lacks random components, which means that with the same data the segmentation will always be exactly the same, which gives total replicability of the results; Experimental data on the improvements obtained in each of these sections are provided below. These improvements consist of:
Innovation: provides the first known solution to perform automatic segmentation of the tear ducts; It is robust, offers C4! results very consistent with the manual method, and 10 exceeds in terms of objectivity, reliability and performance, as set out below.
Sturdiness. The method was developed and tested with the tear ducts of 17 people, normal and with pathology, and once operative the same program was applied to segment another 34 tear ducts of 17 different people, also including normal and obstructed. The program, without the need for any change, segmented the 34 new ducts with a similar performance without any exception and without any incident.
Concordance between the automatic segmentation program and the radiologist. To compare the new automated segmentation with that of an experienced radiologist, 20 tear ducts of people were randomly selected. each of a different person, half with normal ducts and the other half with obstructive pathology. The agreement of the automatic method with the manual has been evaluated in two levels: the coincidence in the pixels considered of the CNL was evaluated by the coefficients of Jaccard, Dice and cosine; the similarity in the area of the section corresponding to the CNL found by each one was evaluated with the Pearson correlation coefficient and linear regression. Fig. 13 shows the value of the overlapping coefficients and 4 ~ 1 value of the Pcarson correlation r = 0.96 (confidence interval [0.91 0.99J, P = 0.05) together with the slope close to The unit of [a regression line certifies that the automatic segmentation provided in the present invention can go to manual segmentation, maintaining the level of quality and gaining in reliability and objectivity.
Performance. The average processing time is about 30 seconds for the segmentation of a channel. This data corresponds to the program that implements an embodiment of the invention in Matlab R2014a environment, running on a laptop with an Intel Core i5 M 540 processor of: 2 cores at 2.53 GHz and with Windows 7 64-bit operating system. Although the program has not been subjected to any optimization pipe and a computer on a basic computer, the execution is significantly faster than that of an experienced radiologist.
Reliability The proposed method is totally deterministic in nature, and lacks heuristic or pseudo-random elements, so the result will invariably be the same.
Objectivity. Being a fully automatic segmentation, subjective factors that affect manual methods are excluded.
Independence and generality This new methodology can be applied directly to any computerized tomography study, in D1eOM or in any other; it is independently of the manufacturer of the Te, of the acquisition technology (multi-cut helical, etc.). of the size of voxel, of the spacing of it! '. cuts and overlaps or jumps between consecutive cuts.
4. MODELING OF THE DUCT (206) The objective of this stage is to automatically generate three-dimensional models of the duct and its axis, starting from the "binary hypermlatrix" B resulting from the segmentation stage.
4.1. 3D modeling of the conduelo (502). Its objective is to obtain various three-dimensional models of the segmented conduit described below.
4.1.1. Direct solid model (504). It is the geometric object formed by the set of "secondary voxels" corresponding to the pixels of all the cuts that have been classified as "conduit" in the segmentation process; the secondary voxel differs from the standard voxel associated with the same pixel in that the secondary voxel may be greater or lesser than the standard, depending on the thickness of the cuts sca mcnor or greater ~ respectfully ~ than the distance between the centers of the cuts ; "secondary voxel" corresponding to a pixel is understood as the solid geometric body defined in the physical space whose center and size of the base in 3D physical coordinates coincide with those of the pixel and whose height is given by the spacing between successive cuts; consequently [the adjacent secondary voxels on any axis do not overlap or leave empty spaces between them.
This procedure to obtain the depth of the advantage advantage over the CI based on "slice thickness" (OO [8,0050), since unlike this it is not affected by the possible overlapping of the cuts or the cxistcncia of space not sampled among them; it also has an advantage over the one based on "Spacing Between SI ices" (00 18,0088) due to the reliability deficit of the latter, given that there is evidence that some manufacturers have introduced erroneous values in that label (Clunie. 2005); for an analogous reason it is also preferable to the one based on "Slic ! Location"
(0020,1041).
4.1.2. Model of the surface that delimits the duct (506). The surface model is obtained from the binary hypermatrix B that contains the sectional masks of the duct by interpolation of the contours of the duct sections in the successive cuts, based on the three-dimensional metric deduced from the tomographic volume metadata , softening the surface resulting from the interpolation, and sclecting the isosupertice as a model of the duct that defines a volume equal to that based on the direct solid modulus. The numerical structure of faces and venices thus obtained represents the surface model of the duct; fig. 16 offers a visualization of it for subject number 17.
4.1.3. Solid fill model (508). It is obtained as the closed set of the points of the space delimited by each surface model.
4.2. Modeling the duct shaft (510). We understand here the axis of the duct as a curve that runs along the duct for its most central part; schematically represents the three-dimensional shape of the duct.
4.2.1. Sequences of fflU (!. "Lareo of the axis of the duct (512). We understand here by the axis of the conduelo the imaginary curve that runs longitudinally through the duct and which is integrated by the central points of the duct. However. Since the duct has been discretized in the tomography with a resolution that in principle is acceptable We will accept the same sampling frequency for the axis We will call the sequence of sampling of the axis to the succession of these points, {CI, c], .... Cp}, being p the number of cuts that intersect the conduit To determine them, two procedures have been used, each of which gives rise to its own sampling sequence 1) Axial sampling sequence obtained as a succession of the centroids expressed in the system of physical coordinates of the axial sections of the duct in successive axial sections: Cr is the point of the physical space that corresponds to the centroid or other central index of the duct mask in the r- th cut intersecting the line (r == 1.2 ..... p); the physical coordinates (xrYr z,) of the c-point, are calculated from the coordinates in the centroid image corresponding to the r-esllno cut, ((, Kr hr), applying the coordinate change matrix of the image space to the physical space obtained from the data contained in
the metadata about the origin and eastern: ion of the fisieo coordinate system and about the spacing of the pixels and the cuts ·: a concrete procedure to obtain said matrix is described and justified in Brett and others (2015). 2) Skeletal sampling sequence, obtained as a succession of the intersection points of the medial planes of the axial slices of the tomography with the axes of symmetry - in the axial direction of TAe - of the voxes belonging to the skeleton :) (pruning) resulting from applying to the direct solid model any of the df ~ thinning or skeletonization methods of images. 3) Sequences of derived sample / inmate. Obtain: nests such as moving average, weighted moving average or other smoothing methods applied to any other sampling sequence.
Each embodiment of the invention can adopt any of the described sampling sequences, which will serve as the basis for the calculation of a multiplicity of models of the axis of the conduit that include: a) the interpolation axes, which are calculated by interpolating the points of the sequence of primary sampling. or from the derived sequence. by any interpolation method; and b) the adjustment axes, which are calculated by adjusting vector functions with IR values) to the points of the sampling sequence, or of the derived sequence, by any adjustment method, the least usual being that of least squares.
4.2.2. Axes of illterpolaciólI (jJ.l). The axis of the duct is obtained by interpolating the sampling sequence by means of a function that will name the axis. Thus the following axes are obtained: - polygonal axis. which is obtained in the following steps: (i) determine the polygonal that joins the points of the sampling sequence of the axis and (ji) prolong each of its extreme segments by an amount equal to half its length. this polygonal being enlarged the polygonal axis of the duct. -the axes if ne, or Bezier. or NURBS (non-uniform rational B-splines), or spline given by the parametric curve obtained by applying the sine or Bezier or NURBS or spline interpolation method to the axis sampling sequence. respectively. -the straight axis, which is given by the rectilinear segment that joins the two ends of the polygonal axis.
4.2.3. Shaft, · adjustment (5/6). The axis of the duct is obtained by adjusting vector functions with values in IR3 to the points of the sampling sequence. or from the derived sequence, by any method of adjustment, being the most common ed of least squares. By way of indication, but not exhaustive, it is worth mentioning the polynomranic axes in which the degree of the polynomial is defined in
I take 3; the parametric expression of the curve that provides a polynomial model of degree
5 1 O 3 is: (xCt) = s Ct) =) yCt) = lzCt) = ao + al z (t) + a2z2 (t) + a3 z3 (t) ba + b, zCt) + b, z'Ct) + b, z3 Ct) t tE (z. - ~ z, + ~) cDl. 1 2 '2 where S is the spacing between cuts and z, and z the coordinates on the z axis of the centers of the initial and final cuts of the channel, the coefficients a, and b, define the axis of each tear canal and its values they are obtained by adjusting l: ada po li nomio by least squares to the respective coordinates of the CI points. e, ..... c "Thus, the calculation of the polynomial axes requires the calculation of the axis of the conduit given by a vector function whose components are polynomials of the lowest degree with which an acceptable adjustment is obtained to the points of the axis sampling sequence.
The axes described above can be calculated both on the direct sampling sequence and on those derived from the la.
15 20 The po lin omic model! I) - (x (t) y (t) z (t)) of the axis of a naso lacrimal canal of subject number 26. expressed in parametric form and in the fisieas coordinate system, is shown in fig. 14 to; Figures 14b and 14c illustrate the method of adjusting said polynomial model in the sagittal and coronal planes; and in both the good adjustment of this function to the empirical points is evident, both for the x coordinates of the centroids of the channel in the axial cuts (~ = 0.9973) and for its y-coordinates (r = 0.9504).
25 30 Fig. 15 shows the values of the coefficient of determination (r) obtained in a sample of 68 nasocriminal channels of 34 people for the polynomial models of grades 2, 3 and 4. The value close to the unit of said coefficient: enough for the Three polynomial models used show the good quality of the adjustment and show that they represent the three-dimensional structure of the channels with confidence and that they also provide a solid basis for the measurements made from them. Thus, the vl: ntajas that contribute these methods are evident. since they provide mathematically precise models and very adjusted to the data; precisely for that reason, they cannot directly compare with the qualitative models described in previous scientific publications or with artistic drawings the shape of the channels.
5. MEASUREMENTS OF THE SIZE OF DRIVING (208)
This stage has the objective of carrying out measurements of various physical dimensions of the three-dimensional duct and its sections based on its segmentation in the tomography and of the correspondence that the metadata of the flamagraphy establishes between the space of the images and the physical space material. 5.1. 3-D measurements (602)
5.1.1. Volumell of the cOllduct based on model only (604). The volume of the duct in the physical coordinate system is obtained as the sum of the volumes of all the "secondary voxels" that form the direct solid model, when the binary hypermatrix B has only three dimensions and all the voxels are the same size , the same result is obtained by summing all the elements of said matrix (which gives the total number of voxels belonging to the conduit) and multiplying the result by the size of a secondary voxel that is given by the metadata and calculated as jodi, ca then.
In the calculation of the volume of a "secondary voxel" the dimensions of its base are given by the DrCOM label "pixel spacing" (0028,0030) and its height by the difference in the third coordinates of "image position (paticlflt)" (0020,0032) in two consecutive cuts; This calculation procedure must be adapted with geometric criteria in the event that the coordinate system indicated by the director cosines on the "image orientation patient" label (0020.0037) is different from the canonical one.
The effect of partial volume on the periphery of the duct is a known source of error in the evaluation of the volume of a tomographed duct. His anisotropic nature. together with the nonlinearity with respect to the proportion of each substance in the voxel. complicates the solution. Here we have opted for statistical cancellation, based on the central limit theorem (Feller, 1971) And justified because, although the number of voxels on the outer face of the channel is small in a cut, however in full volume It is sufficiently numerous to justify the applicability of this theorem.
When the duct is open at one or both ends. as in the case of tear ducts, there is another source of error to be corrected, which is linked to the ends of the duct and which we call "wedge-shaped", It consists in that it is usually considered that the duct begins at the first cut in which the Camama of its section forms a sharp and tennin curve in the last cut that meets the same condition. The error occurs at each end when the plane of the lomography is not perpendicular to the duct and manifests itself in that the wedge of the duct between the plane orthogonal to the duct that actually delimits it at that end and the plane of the tomographic cut is erroneously included or excluded in the volume calculated by the proposed method. The error derived from the wedge effect (which may be due to excess or default) is not compensated for errors of the opposite sign and therefore cannot be considered corrected by applying the central limit theorem. To correct it, it is essential to add a correction to the volume obtained before, which is obtained in this way: a) determine the outer plane of the extreme cuts that contain the duct and that appear in the general parameters; b) determine the plane most distal and orthogonal to the axis of the duct whose intersection with the surface that delimits the duct is a closed curve; c) calculate the volume of the wedge between said outer plane, said more distal plane and the outer surface of the duct; d) subtract said correction to the volume of the rest of the duct if said wedge was already contained therein
or add it in an ordinary case: and d) if the duct is also tubular in shape and curved in the extreme air, in lances repeat all the previous steps also for that other end.
Advantages. The method proposed here for the calculation of the volume has the following advantages over the existing methods: a) it is supported by an automatic and objective segmentation, instead of manual and subjective: b) takes into account and corrects the errors due to the "wedge effect ", effect ignored in previous scientific publications; and c) correct errors due to the effect of partial volume more easily than others such as those based on sub-pixel precision (Huertas and Medioni, 1986) or on correction techniques of the peTimctro and the ark such as those proposed by Pran (1991, chap. 8.2).
5.1.2. V () lumen of the conduct () surface model and fill () (606). The volume of the duct delimited by the surface model or the solid filling model is obtained by using numerical integration techniques or finite elements.
For the two exposed methods for calculating the volume of the duct, and when the duct is open at one or both ends, a correction must also be applied to the volume obtained so that for each end it is silenced so a) determine the outer plane of the extreme axial cuts that contain the duct; b) determine the plane most distal and orthogonal to the axis of the duct whose intersection with the Sup (: area that ends the duct is a closed curve: c) calculate the volume of the wedge comprised cntrc the aforementioned outer plane. the cited plane more
distal and outer surface of the duct; d) subtract said correction to the volume of the rest of the duct if said wedge was already contained in it or add it otherwise.
5.1.3. LOlJgilu (¡del cOlJduclo (608). The method to obtain the length of the duct depends on the type of axis taken as a model: a.xial, straight, polygonal and poxomic of degree n. The length of the duct is obtained by calculating the length of the curve that represents the axis adopted as a model and between the planes orthogonal to the conduit that delimit it at its two ends.
Figure 18 a) shows the means of the measurements of the length of 68 tear ducts. In what we know are the first measures obtained in a fully automated way. 5.2. 2-D measurements (6/0)
5.2.1. Area! Ieccioml¡ axirll (612). In a tomographic image the axial sectional area of the conduit in the physical coordinate system is obtained by adding the areas of the pixels marked as "conduit" in the corresponding binary mask of the conduit, the area of a pixel being equal to the product of its based on its height expressed in physical coordinates by tomography metadata. When all the pixels are of the same size, the same result is obtained by adding all the elements of the binary matrix corresponding to that image in the binary hypernatrix B (which gives the total number of pixels belonging to the conduit in the aforementioned tomognitic image) and multiplying the result by the surface of a pixel expressed in physical units and deduced from the metadata; When the image is generated in the D1COM standard, the surface of a pixel is obtained by multiplying the base by its height expressed in the D1COM tag "pixel spacing" (0028.00: 30).
To obtain the axial sectional area of the duct at any point on the axial axis, the image of the axial cone at that arbitrary point 2fJ of the z-axis of the tomography must first be obtained, together with its metadata. and the corresponding binary mask of the duct section in said image. To do this, follow these steps: a} If the value of 20 is equal to the z-code of any of the cuts of the original tomography, then the desired result will be formed by the image. the metadata and the binari mask: a corresponding to said cut; b) in another case, multiplanar reconstruction must be calculated keeping the orientation and spacing of pixels and cuts the same, but matching Z (.I the z coordinate of one of the eons: apply to multiplanar reconstruction the preprocessing method and segmentation; the desired result is given by the image of the zr¡ cut resulting from the multiplanar reconstruction, along with its metadata, and the corresponding binary mask in the binary hypermatrix resulting from the segmentation. Calculate the sectional area of the duct in the cut that goes through Z / J following the same procedure as for the cuts of the original tomography.
The axial area has been used profusely, c, but it has the problem that it grows artificially as the cutting plane moves away from the position perpendicular to the axis of the channel. In addition, the errors induced by this reason are systematically positive. The inevitable consequence of this is that the area of the axial section of the duct may not be indicative of its actual section. Having the axis of the duct obtained by the innovative procedures described above, now allows solving the problem posed by the axial area and assessing the actual sectional area by estimation and even by direct evaluation.
5.2.2. Section section! Orthogonal estimate (61./). It is obtained by multiplying the corresponding area in the axial section by the cosine of the angle between the axis of the tomography and that of the duct; that is to say. the estimation of the orthogonal sectional area from the axial is obtained by projecting the axial area on the plane perpendicular to the e; e: therefore in the section k the sectional area orthogonal to the channel At: is given by A ~ = AkCOS PL where Ak it is the axial area of the cut kyp ~ it is the angle that forms the axis of the channel with the perpendicular to the tomographic cutting plane.
To calculate the orthogonal sectional area. Estimated in the plane perpendicular to the axis at any point p thereof. the process is conceptually analogous although it requires more steps, which are the following: a) calculate the tangent in p of the mathematical function adopted to model the axis; b) calculate the multiplanar reconstruction maintaining the orientation and structure of the original cuts but then loosening them on the axis: so that the cut closest to point p happens to contain that point; c) obtain the binary duct mask in the image of the cut that passes through p by applying to the multiplanar reconstruction the methods of preprocessing, segmentation and calculation of the axial sectional area for the cut that passes through p; and e) calculate the estimated sectional area by multiplying the aforementioned axial sectional area by the cosine of the angle between the axial axis of the tomography and the tangent obtained in the first step.
5.2.3. Orlogollal sectional area (616). The method to obtain the real orthogonal sectional area (not estimated) in the plane perpendicular to the axis of the duct by an arbitrary point p, consists of two phases: the first is to obtain the binary mask of the duct section in the orthogonal section.
to its axis at point p; and the second calculate in said binary mask the axial sectional area. thatby construction it coincides with the sectional area orthogonal to the axis of the duct.Jj obtaining the image of the orthogonal cut to the axis of the duct at an arbitrary point p ofthis, together with the corresponding metalclales and the corresponding orthogonal binary mask.
5 includes the following steps: a) calculate the tangent in pa the mathematical function adopted as the axis model: b) determine the equation of the plane perpendicular to the axis of the duct in p: c) calculate the multiplanar reconstruction containing said plane perpendicular to the axis in p, juma with the corresponding metadata. and obtain the orthogonal sectional image of the duct at point p along with its metadata. what will it be equal to the image of the cut corresponding to said
10 point in said multiplanar reconstruction, along with its metadata; d) apply to the multiplanar reconstruction the segmentation method of the invention and obtain the orthogonal binary duct mask t! n the cut pt = rpendieular to its axis at point p that will be equal to the binary matrix corresponding to point p included in the binary hyperpentry matrix, 3) The orthogonal sectional area is obtained by adding the pixel surfaces that the
15 corresponding orthogonal binary mask: it is indicated as "conduit" in the selected coordinate system: the surface of a pixel in the physical coordinate system is given by the product of its dimensions indicated in the metadata. which in DICOM gives the label "pixel spacing" (0028,0030); in the image coordinate system the [size is arbitrary and is usually taken as the unit in both axes.
Advantages. This method improves the one proposed by Ramey et al. (2013) since it is automatic instead of manual; it relies on objective data on the orientation of the channel instead of subjectively estimating it; it is faster, and consequently it can be evaluated in a high number of points. situation in which the manual method results in the practice little viable.
25 l a Fig. 18 b) collect the means of the measurements of the axial and orthogonal sectional areas estimated based on the polygonal axes. and the polynomials of 21. JO AND 4th grade. All of them have been obtained without any manual intervention, 10 which represents a novelty in what we know.
30 The minimum sectional area, axial or orthogonal, represents the minimum value that said area reaches along the sampled points of the channel. The depth of the minimum, axial or normal area. is the length of the axis measured from the beginning of the channel in the orbit to the point where the corresponding area reaches its minimum value.
35 5.2.4. Axial orientation J 'orlogon (ll of 1'0 section (618).
The orientation of the axial section of the duct given by the orientation of its maximum diameterin the chosen coordinate system.The orientation of the orthogonal section to the axis of the duct may be different and is given bythe orientation. in the chosen coordinate system. of the maximum diameter.
5.3. Measurements 1-0 r620;
53.J. I. Length of diameter max; mIJ or ml "imo in the axial section, (622) calculated as the maximum or minimum distance - respectively - between the outer edge of all the pairs that can form the pixels located on the periphery of the section duct
5.3.1. Maximum and minimum orthogonal diameter length in orthogonal corle f62.J)
calculated in the coordinate system chosen as the maximum or minimum distance respectively-between the outer edges of all the pairs that can be formed by the pixels located on the periphery of the duct section.
6. MORFOMETRtA
This stage has the objective of obtaining descriptive indexes of the canal's duct and its sections. from the information provided by the numerical models of the same obtained in an earlier stage (fig. 7).
6.1. Axial section morphometry. "And orthogonale.'i (702).
6.1.1. Moments (70-l). The moments (! Standards (mfll), mlfl. Mili, m il. M} (J, mllJ. MJ11. MnJ. Mi}. Mu), central (~ I. J.111. ~ LO. J, lol. ~ .3IJ, J, loJ. ~ ¡Il. J.1} /), standardized central units (nno, n //. N 2D, n D). N J (J, nllJ. N I], n, l) , and invariants (MI. Ml. MJ. M.¡, M5. Mr .. M7) are calculated in the coordinate system of the image and in the physical from the contours of the segmented sections of the binary hypermatrix, uti Using the formulas provided by Levine (1985), Gonzalez and Woods (2002) and Hu (1962), all were calculated in physical coordinates.
6.1.2. Eccentricity (706). It represents the degree to which the section of the channel moves away from the circular shape: when the section is a cyn: perfect ulo. the eccentricity will be worth O: as the section takes more elongated forms the eccentricity will take positive values increasingly closer to . The axial or orthogonal eccentricity of a flat figure is calculated as the ratio between the distance between the foci and the length of the major axis of an ellipse with moments
non-centralized second-order exchanges equivalent to those of the duct section in the
image of the orthogonal cut to the axis of the duct, or by any other available procedure. The axial eccentricity of the orthogonal is differentiated, as evaluated in an axial section or one orthogonal to the axis of the duct, 6./.3. Elongation (708). It has an interpretation similar to eccentricity. The axial and orthogonal can be distinguished. according to YES! evaluate axially or orthogonally to the axis of the duct. In both cases, the difference between the unit and the ratio between the minor and major sides of the enveloping rectangle of the axial or orthogonal section of the duct, or by any other available method is applied by applying the metric described for this in the metadata.
6.1.4. CirclIlarity (710). It can be understood as the opposite dimension of eccentricity and elongation. The proposed formula Zunic, Hirota and Rosin (2010) is used to calculate it
where C (S) represents the circularity of the region S and f.1¡k represent the central moments of order} and k. 6./.5. Inrlice.'i tle mor / ome / ría geometr; ca 2D (712). The method for carrying out geometric morphometry amilysis of the fonna of any section of the conduelo consists of a) calculating from each binary sectional mask a) (lal or orthogonal the corresponding line of axial or orthogonal sectional conrorno -respectively- in the chosen coordinate system, a line that as such has a null thickness and is integrated calculated as the set of contact points with the background of the pixels that belong to the conduit; b) generate the marks ("Iandmark.o) ¡") And the matrices of reference points in physical coordinates established in the metadatas of the contour of the binary image of the conduit in the orthogonal cut to its axis, and represented in the corresponding binary sectional image; and e) calculate one the plurality of indexes of the orthogonal section that offers the geometric morphomctrics calculated directly on the line of conlnrnn of the conduit OR starting from marks or "Iandmarks" derived from it. in its same coordinate system. and that will include analysis of Procrustes (Bookstcin. 1991: 1997). analysis of allometry (Klingenberg, 1996), and analysis of symmetries and asymmetries (Klingenberg and Mclntyre. 1998).
Finally. the obtaining of statistics of each morphometric variable of the axial or orthogonal sections (depending on the type of cut) by: a) the repeated measurement of said variable at all points p, of the sampling sequence of the axis, or another set of points of the axis that is of interest: b) the calculation of the statistics of the measurements obtained. which include one or more of the following indices: i) the indices of central tendency and variability, ii) their extreme values. minimum and maximum; iii) determine the axis points where
these extreme values are found; and iv) obtain the curves of each variable that express the measured value in ordinates and in abscissa the point of the axis of the duct or the axis of the tomography. 6.2. MorComelria duct (714)
5 6.2.1. Specific indexes on the axis of the duct (7/6) This section provides methods for calculating morphometric indices at any point on a previously calculated axis and that include b) the orientation at each point of the axis and the average orientation at the points of the sampling sequence; c) instant curvatures. average and maximum; d) instantaneous, average and maximum torsions; e) the location of the maximum points
10 and minimum curvature and torsion; l) Euclidean distances from said singular points to each end of the axis; and g) the straight distances and along the axis between said singular points and every ~ xtr ~ mo dt! 1 t: j ~. to Orientation of the axis of the duct at an arbitrary point. When the axis is a function, differentiable the orientation of the axis at any point p is given by its tangent Have said
15 point: T (p) = (r '(p) 1] 1 r' (P) 11; when it is not, as with the polygonal axis at the points of the sampling sequence of the axis, then it is calculated taking as free vectors the segment that arrives and the one that leaves the vertex and taking as an orientation of the axis the one of the vector sum of both b) Curvatu ra of the axis of the conduelo, in an arbitrary point. For polynomial axes,
20 dependent on parameter I then the curvature c is obtained using its first and second derivatives, r 'and r ", with the formula (/ _ 11.' (1) X," (1) 11
) -III "(tll ~ For the polygonal axis, which is not differentiable at the points of the sampling sequence, segment l ', which reaches the k-th vertex and the one leaving it, are taken as free vectors Sl + /, and from
25 of them are calculated as the index of the curvature at those points, the angle between both vectors has been obtained using the formula u ~ (; us · 1 <Stoi "t + ¡>
15. o: IIsl <+ II e) Shaft torque. Torsion reflects. to what extent the major axis of the channel rotates as it goes from one end of the duct to the other. When the axis is differentiable as in the case of
30 polynomial torsion t is calculated using the formula. '(1) · (1 "' (1) x 1" "(1))
(
Tt)
- 111 "(1) x 1" '(1) 11'
where, ',, "and,'" are the first three derivatives of the axis; when it is not differentiable then it
Calculate by numerical methods.
6.2.2. Scitional and punctual indices ~ the long axis (718). They are obtained by: a) repeated performance of the measurements of said variable at all points p. of the sequence of sampling of the axis, or other set of points of the axis that is of interest; b) the calculation of the statistics of the measurements obtained. which include one or more of the following indices: i) the indices of central tendency and variability, ii) their extreme values. minimum and maximum; iii) determine the points of the axis where said extreme values are located; and iv) obtain the curves of each variable that express the measured value in ordinates and cn abscissa the point of the axis of the
10 duct or axis of the tonography.
For each variable that has been measured at several points along the axis, such as the
sectional area or curvature. the curvn quc is expressed expressing the measured value as a function of the point
of the axis of the duct or of the tomography in which the measurement was made. For example. for him
In the sectional area, the value of said variable in each cone will be placed in first order and in abscissa the value of the depth coordinate z of each cut (in mm); and secondly, the value of said variable in each cut and in abscissa the length of the section of the axis that goes from its origin to the point from which the measurement comes is placed.
20 For each sc curve, its mcdiantc shape will be adjusted for adjustment of polynomial or other curves and the equation of the curve, its characteristics and the interpretation of these in tenninos of the duct fan will be recorded; for example, if the curve corresponding to the orthogonal sectional ark index is linear increasing or decreasing, it will indicate that the duct is shaped like a cone with its wide part towards the end of the axis or towards the beginning. respectively; If it is constant. indicate to
25 that the duct lacks duct narrowing (uniform section); If it is parabolic. It will indicate that the duct has a double funnel shape (with a widening in the central part) or an hourglass shape, depending on whether the sign of the second degree term coefficient is positive or negative (quote).
30 6.2.3. Location on the axis of the point. S sylgular (720). The maximum and minimum values of each variable measured at several points along the axis allow defining the point of the axis where it is located as a singular point. that characterizes the duct and that can be identified by its 3D coordinates or by simpler indexes such as its distance to one of the ends of the axis. Thus it is possible to define new variables that characterize the structure of the duct indicating
35 where are its singular points. such as the depth at which a maximum or minimum sectional area, axial or orthogonal, is located; or the depth at which the axis maximum curvature point is.
b) the orientation at each point of the axis and the average orientation in the pumos of the sampling sequence; Instant curvatures. Medium and maximum: instantaneous torsions. average and maximum; the location of the points of maximum and minimum curvature and torsion; ) the Euclidean distances from said singular points to each end of the axis: and) the straight distances and along the axis between said singular points and each end of the axis.
In general, for all previously defined indexes that are measurable at different points on the axis, their statistics are obtained by means of: :: a) the repeated realization of the measurements of said index at all points p, of the sampling sequence of the axis, or other set of points of the axis that is of interest: b) the calculation of the statistics of the measurements obtained, which include one or more of the following indices: i) the indices of central tendency and variability. ii) its extreme values. minimum and maximum; iii) determine the points of the axis where these extreme values are found: and iv) obtain the curves of each variable that express the measured value in ordinates and in abscissa the axis of the lomography or the rectified axis of the axis of the duct.
6.2.4. 3D geometric morphometry indexes (722). The method for carrying out geometric morphometry analysis of the three-dimensional shape of the duct includes the methods to: a) obtain the landmarks and reference point lines in the physical coordinate system established in the metadata, from the 3-D model of sectional contours or the surface model or the axis model of the duct: and b) calculate a plurality of 3D geometric morphometry indexes based on said marks and reference matrices. which include Procrustes analysis (Bookstein. 1991: 1997). allometry analysis (Klingenberg, 1996), symmetry and asymmetry analysis (Klingenberg and Mcl ntyre, 1998).
6.3. Dynamic or evolutionary study (7.24). It consists of the repeated application of the previous analyzes to successive tomographies of the same conduit obtained over time, which also includes the methods to a) calculate evolutionary indices of each variable and b) generate the curves and other tannin models representative of the evolution of these variables throughout
weather.
6.4. Creation and study of the lomo ~ r'afia rectificltda, following these steps: a) obtain the image of the orthogonal section to the axis of the product at all points of the sampling sequence of the axis; c) position each image thus obtained in a plane perpendicular to the axial axis that cuts it by the z-coordinate corresponding to that image and moving it in the indicated plane so that the cenlroidc of the duct section co-enters with the axis, obtaining thus the rectified lomography of the duct; d) record it in permanent memory: e) represent it graphically for visual and morphometric analysis.
7. RESULTS (21 2)Once the automated analysis is finished, all the results are permanently recorded.obtained and made accessible to the user in local or rem mode using the mostappropriate in each case described below:
1) For three-dimensional models of axes and ducts, the method includes a) permanently recording the models obtained: b) making them accessible to the user. visualizing the models in virtual reality devices, screens or paper, in their entirety or by cuts; isolated or integrated in a volumetric rendering or other representation of its environment; In video mode. static figures graphics. or interactive figures allowing the user to select the observation point. the colors and nature and position of the light sources; c) print in 2D or 3D of the results chosen by the user; and d) generate digital copies on fixed or removable storage media.
2) for 2D graphics and numerical results, the method includes a) recording the obtained graphical and numerical results in permanent support; b) make them accessible to the user, displaying them on screens. paper or other volatile or permanent physical media. in its entirety or in parts: in video mode. static figures, graphics, or interactive figures; c) print the results chosen by the user; and d) generate digital copies on fixed or removable storage media.
Advantages. The method proposed here for (; alculating the volume of the tear canal from its segmentation in sections a) (iales has the following advantages over that proposed by Estes et al. (2015) a) is based on automatic segmentation and objective, instead of manual and subjective, b) it makes a statistical cancellation of the segmentation values 20 instead of accumulating the subjective biases that manual or assisted segmentation has been able to introduce): e) corrects errors due to partial volume effect: d) correct the enores due to the wedge effect. that in previous publications is ignored; and e) generates measurements with an accuracy of fractions of mm), instead of whole numbers of mm3 • If compared with previously published methods, the main advantages are: a) it is realist, since it evaluates the actual volume of the channel according to its rcal form, while the previous methods are limited to estimate its volume based on approximations of its shape given by regular geometric figures; b) it provides greater precision in the measurements: e) it is an objective method, while the others are subject to the subjectivity of the operator who decides in each case at what points and between which points he takes the measurements: and e) is faster .
All the analyzes described have been carried out independently for each of the 68 lacrimal channels studied and in less than 10 hours of calculation in a laptop. In what we know. The method of this invention has allowed for the first time. thanks to the automation of segmentation, first define an innovative methodology for the study of tear ducts, and secondly obtain massively objective measurements in a number of variables close to 2010 with which new ways for three-dimensional study are opened of structures such as tear ducts. The method thus represents an innovative tool to promote the advancement of knowledge in the areas of ophthalmology and plastic and constructive surgery. as well as in other areas that until now posed similar problems. It therefore provides a clear advantage over previously existing methods, focused on few variables, and often based on approximations and always with measurements whose reliability was reduced by the operator's subjectivity when performing them manually.
The specific advantages of each of the methods presented over the existing ones have already been indicated in their corresponding sections. In global terms and in relation to the advantages of this invention. It should be noted that thanks to its automation a) it allows to obtain modeling and measurements with rapidel ~. precision and objectivity: b) is applicable to any tomogenic modality (CT, NMR, PET, SPECT, OCT, focal microscopy and excitation fluorescence with one and two photons ...); c) is independent of the manufacturer of the tomography equipment and its technology (helical, multislice, confocal, fluorescence, etc.); d) whatever the procedure followed to generate the tomognific volumes, the IOmognifica reconstruction core, the thickness and separation of the cuts and the size of the pixels ~ e) can be used in the case of overlapping cuts or with spaces not scanned between them and even with sequences of variable thickness and spacing cuts: 1) eliminate the random variability that manual methods introduce; g) do not require any assumptions about the shape of the channel section in each case; h) take more advantage of all the
information provided by the spelling, lament in the image data as in their headers: and are applicable in medicine. and other scientific and industrial areas.
DEVICE
A storage device penlHmente program readable by a computer. which tangibly incorporates a program of executable instructions in a computer. a digital signal processor, an application-specific integrated circuit, a microprocessor, a microcontroller or any other form of programmable hardware capable of performing the steps of the method described in this invention.
Figure 1 offers a diagram of the structure of modules that integrate the device of an embodiment of the invention consisting of a computer (100) having its own calculation unit
or CPU (122). its working memory (124), its input devices including mouse and keyboard (126), its output devices including screens with the ability to display data and images (128), a non-portable file system for storing data and Resulted (130), a non-volatile and tangible medium for archiving d <: programs (110) that can be read and executed by the CPU with its working memory; said non-volatile medium. tangible and readable by computer for program archiving (110) contains the program modules in which the methods of obtaining the tomography are materialized (11.1). preprocessing (112), segmentation (113), modeling (114). measurement (115), morphometry (116) and recording, presentation and use of results (117); all the modules that make up the computer are connected to an internal communication bus (132) that puts them in two-way communication with each other and with the CPU and working memory; also through the external communications card (134) the computer also communicates with the external communication networks (136), thus being able to obtain directly through dt: they the t: omographs from the tomographic equipment
(108) Or of the PACS OR files where they are stored (106) and even managed by users who interact from their own devices through the Internet (104) or from a remote workstation (102).
SYSTEM A system comprising the device described above and a web server. connected by a network to said device and another to intemet or other telematic networks; said server a) pennile Connect to remote users, receive the sending by the remote user of a request for analysis and modeling and the tomographic volume to be analyzed, together with the complementary infoOllation; b) transfers the request and the tomogenic volume to the aforementioned device to perform the analysis of the three-dimensional and automated dimensions and morphometry of the duct; e) receives the results in electronic format from the device; d) transfers them to the applicant through the lelemmatic or other channel that he has selected: and e) informs of the completion of the process to the provider's information systems.
BRIEF DESCRIPTION OF THE MBUJOS
To complement the description that is being made and in order to help a better understanding of the characteristics of an embodiment of the invention, a set of drawings is attached as an integral part of said description, where it is illustrative and not limiting, It has been: prest: ntado lu next: nte:
Figure 1.- Shows a diagram of the structure of modules that integrate the device of an embodiment of the invention, as well as the functionality of each of the modules.
Figure 2.- Shows the steps that integrate the method of an embodiment of the invention.
Figure 3. · Shows the main stages of the preprocess (202). which includes the necessary procedures to create and record all the data necessary to carry out the automated segmentation of the duct into the tomography.
Figure 4. · Shows the steps included in: 1 automated segmentation procedure (204) of an image. as well as the steps that each stage includes until first obtaining the segmentation of the conduit in that image; and finally obtain the three-dimensional segmentation of the duct by applying the same procedure described for an image to all the cuts of which the tomography consists.
Figure 5. · Shows the main stages of the modeling procedure (206) of the newly segmented three-dimensional duct. and that includes methods for modulating the segmented conduit (502) and its axis (510).
Figure 6. · It shows the main stages of the measurement procedure (208), including the measurements of the three-dimensional duct (602), its sections (610) and its axis (620), as well as the stages that consist each of the above.
Figure 7.- It shows the main stages of the procedure to obtain indexes of the duct shape (210). and which include morphometric indexes of the sections (702), of the three-dimensional conduel (714) and of their evolution over time (724).
Figure 8.-Image of a tomogenic cut of the head: the arrows at the top indicate the location of the sections of the two tear ducts. which are the object to be segmented in the embodiment of the invention.
Figure 9.- Shows the preprocessing stage consisting of the delimitation of the region of interest in the image (306): a) partial and enlarged view of the part of the image in which the section of the right tear duct appears: b ) region of interest indicated by the thick black line, consisting of a polygonal of few vertices that encompasses not only the nasolacrimal canal but also a strip of the surrounding substance. in order that the delimitation of this region of interest does not condition the automatic segmentation process; the other fine lines introduced in the image represent line: reference ils to speed up and verify the layout of the region of interest; c) mask of the region of interest delimited: d) area of the image marked by the mask of the region of interest in which the automatic segmentation process will focus.
Figure 10.- Shows for the segmentation stage (204) the result of the feature extraction stage in the region of interest (414): a) feature extracted by the bi-threshold operator: in the white pixels of the part central the trait takes the value 1 (suggesting "channel"), in the gray pixels that mostly surround the whites the trait is worth 0'5 (reflecting "doubt") and in the dark gray pixels the trait takes the value Or (suggesting "bone") b) feature extracted by the centripetal directional derivative operator: in the lightest pixels the slope of the derivative is greater and positive, in the darkest the slope e: s large and negative and in gray next to zero: and c) trait extracted by the Laplandian operator: in the darkest pixels in the Laplandian it is negative and large, in the lighter ones the other way around. and in gray near zero.
Figure 11.-Shows for the segmentation stage (204) a) the provisional duct mask, resulting from applying the classifier to the pixels of the region of interest (416). b) the rectified mask resulting from the rectification step (418), c) the final binary mask resulting after the iterations (422), and the region delimited by said mask in the initial image.
Figure 12.-Correlation between an expert radiologist and the program described in this invention in terms of the number of pixels that form the section of the nasolacrimal canal in a random sample of 20 cuts belonging to head scans of as many people.
Figure 13.-Statistics and indexes of coincidence between the segmentation performed manually by an expert radiologist and automatically performed with the program described in the preferred embodiment of this invention.
Figure 14.- Third degree polynomial model for the right channel of subject number 26 of the sample: a) parametric equations of the model, b) adjustment of the model to the points in sagittal projection and e) in coronal projection; each point represents the centroid of the section of the lacrimal naso canal found in each section by the program described in the preferred embodiment of this invention; Its coordinates in mm are referred to the coordinate system of the physical space.
Figure 15.- Adjustment indices of polynomial modelDs of grades 2 ' 3 ° Y4 ° to the centroid wings of the sections of the canal: data from a sample of 68 naso-tear channels of 34 people.
Figure 16.- Shows two views of the surface model of 105 two naso-tear ducts for the
subject number 20 of the sample: in the bottom view a curve has also been included in each channeldc dark stroke representing the model of the axis of the channel.Figure 17.-Elements estruclurantcs in rorma of L.Figure 18.-Averages of lengths (a) and sectional areas (b) obtained automatically withthe software of the preferred embodiment of this invention in 68 tear ducts of 34 subjects.
PREFERRED EMBODIMENT OF THE INVENTION
While the invention is described and illustrated in a preferred embodiment. specifically in its application to the segmentation and three-dimensional analysis of the tear ducts in computed tomography for clinical use, the invention can be applied and produced with many different configurations. It is represented in the drawings, and will be described herein in detail. A preferred embodiment of the invention, well given that the present description is to be considered as an exemplification of the principles of the invention and the associated functional specifications for. His construction; It is not to be understood that the invention is: limited to the illustrated embodiment but also covers equivalent embodiments, even applied in other areas such as the study of the mandibular canal
or in industrial applications. Those skilled in the art will conceive many other possible variations within the scope of the present invention.
The preferred embodiment described in this section is to quantify in an automated way the metric and morphological characteristics of [the tear ducts of living persons in head tomography images obtained with tomographic equipment and clinical use protocols.
1. OBTAINING THE TOMOGRAPHY (200)
1.1. Acquisition of the lomographic volume. The tomographs were obtained from a PACS in D1COM format. They were generated with a team for clinical use and without contrast, with the subject in supine position and the image plane perpendicular to the table; helical acquisition; standard or bone reconstruction core; Jada 0.3 to 0.4 mm square pixels: spacing between cuts 0.6 to I J mm.
1.2. Maximize the number of cuts that include the channel. In the lomographs used, the angle between the axial cutting plane and the tear ducts is not too far from [erpendicular. So at this stage the transformation applied was identity.
2. PREPROCESS OF THE HEAD TAC (202)
2.1. Create the files with the image data and metadata (302). Image data and metadata were separated using hoc programs; to the image data the transformation H = a + bT was applied with the values of a and b that for each image provide the DICOM tags (0028, 052) Y (0028,1053) respectively in each TAC; In this case, all indicated a = -1024 and b = 1; with them, the hyperpeller V of 51 2x512xt was created, 1 being the number of cuts in the volume; and with the metadata the matrix M of tx2 was created. putting in the first column the text indicating the meaning of the values placed in the second column of the same row.
2.2. Normalization of the fi ~~ ura-bottom gradient (304). In this embodiment the figure-bottom gradient is already notualized in the original tomography.
2.3. Define the region of interest (ROl) of each channel in sludge 105 cuts (306). For the delimitation of the region of interest in each image, a computer-assisted delimitation was chosen to avoid inadvertent errors. To do this, the image of the channels (x8) was enlarged.
smoothing it through a low-pass Gaussian filtering to eliminate the visual mosaic effect, adjusting the window level to 300 Hounsfield (UH) units as Maatman proposes
(1986, p. 15) And its amplitude at 1000 UH to make the contours of the canal more visible (fig. 9 a); to guide the path, the laplacian zero crossing lines and the isocontorns with UH = 176 were superimposed on the image in red and green (both in gray curved lines in Fig. 9b).
With these aids, the ROl was manually delimited by means of a polygonal with sections of variable length at the discretion of the operator, always by the bone area, but slightly away from the apparent contour of the channel: the prohibition of crossing the red and green lines is also established when they surround it, in order to ensure that the manual delimitation of In RDI will not interfere with In automatic channel delimitation: a polygonal thus obtained is represented in fig. 9b along the thick black line. Next, the binary mask of the region of interest is generated (fig. 9c) that defines the ROl and that limits the area of the image to which the duct segmentation process is to be restricted (fig. 9d).
2.4. Determine general parameter values (308). The following general parameters were included: a) the number of the first and last image to segment; b) the criterion for stopping the segmentation iterations with each image, based on the sum of quadratic differences between the current mask and the previous ones and requiring at least two iterations in each image and establishing termination when the resulting segmentation mask is equal to aira obtained previously, c) [os parameters of the multiumbral operator: Bm ", = 600 whose value represents the minimum anticipated value for pixels not belonging to the channel; C" ,, "'' '50 and Cmu = I 00, which respectively represent the lower and upper dimensions for the range of the maximum expected value in all the cuts for the pixels belonging to the "channel" (the values of Cma ~, Cm, "and Bmon have been established by an expert based on the available data in the
specialized scientific literature); d) the parameters common to the operators derived from discretionary and discretionary derivative: pixel jump equal to l and vertical and horizontal distances between pixels given by the D1COM tag "pixel spacing" (0028,0030); e) the shape and size of the structural elements for the mathematical morphology operations detailed in a later section; and f) the addresses of the input and output files.
3. CHANNEL SEGMENTATION (204)To obtain the segmentation of the channel in each cut and in the volume, the described method is appliedto contil1uación and that includes the following stages (fig. 4):
a) Evaluate and record the values of the specific parameters of the image necessary for the next iteration (412), and that induce a) determine the region dI! reference, that for the first iteration is the RO l, and for the sf: gunda and following is determined by the mask of the channel obtained in the previous iteration: b) determine the area of the region that in the physical space 5 corresponds to the reference region of the image. using the information on the size of the pixels contained in the metadata from the orCOM tag "pixel spacing" (0028,0030); e) determine the physical coordinates (e, el) of the central reference point. which is the ccntroid of the reference region (RR) corresponding to each iteration; d) determine the size of the circularly squeezing elements used by the
10 mathematical morphology techniques - erosion and dilation - depending on the size of the reference region (RR): the radius (measured in pixel: s) is worth 1 if area (RR) <300, 2 if 300. $ ilrea (RR ) <600. 3 if 600 :: Sarea (RR) <900, and 4 if 600 $ area (RR).
b) Calculate the associated feature neighbor 11 each pixel of the ROl (414). This stage consists
15 in the application of a filtering for noise reduction and of the 3 operators extracting selected features, each of which obtains the corresponding coordinate in the vector of features associated with each pixel of the image. using the general and specific parameters. The stage includes:
20 1) A Gaussian low-pass filtering to reduce noise in the tomography images. 2) Apply the biumbral operator whose input is the resulting image of the Gaussian filter, and which calculates the first coordinate of the feature vector for all ROL pixels. The threshold uses two thresholds. l / y t2 calculated like this: 11 = mjn [100, max (50. Mrl)] and t] = (600+ 1,) / 2. where M rl is the median of the UH values in a 3x3 submatrix of the image
25 centered on (c, c¡): the value assigned to the first feature and saved for each pixel in the first coordinate of its feature vector. it is I (interpretable as "channel") when X,} <11 • O ("background") when xy> 1] Y 0.5 ("doubt") otherwise. CIII4f and C ", m are the general parameters that delimit the expected range of values for the minimum value in the pixels that belong to the channel; and 8 mm is another general parameter whose value represents the minimum value that a pixel that cannot take can take.
30 belong to the channel. Fig. lOa graphically illustrates in each pixel the values assigned in a cortc: the lightest region in the central part represents the pixels whose feature vector has received in its first coordinate the value 1, the light gray region that mostly surrounds the previous one corresponds to the pixels whose first coordinate took the value 0.5, and the dark gray region represents those that: hiln received the value O; the rest of the pixelcs located in the black zone have not been analyzed by the luminaire and have directly received the value O for being outside the ROL 3) The differential operator used to obtain the value of the second feature, and that will be stored in the second coordinate of the ve: ctor of features. he is the discrete Laplacian operator; it is calculated in the physical coordinates and with adjacent pixels, and the result obtained gives for each pixel of the ROl the value of the second coordinate of its feature vector. The implementation of the V2G • I convolution calculations has been carried out in the frequency domain by
the product of the Fourier transforms of V ~ G and the image I, in order to achieve high filtration accuracy without having to resort to large convolution cores that lead to excessive calculation times. The value thus obtained in each pixel constitutes the second coordinate of the feature vector. The lig. IOb illustrates the value of the Laplacian in each pixel: the clearer, the greater and more positive Laplacian; The darker Laplacian bigger and negative.
4) The value of the third trait calculated is given by the derivative of the UH values in the direction of the central reference point. which in this embodiment is the centroid of the reference region in each iteration. Fig. The graphically shows the value obtained on each pixel of the ROl: the brightness indicates the value of the centripetal directional derivative: in the light zones the derivative is positive, and in the dark ones negative.
Both the discrete derivative and the Laplacian are based on contiguous pixels and use the distance between them indicated for the two axes by the metadata from the OICOM tag "pixel spacing" (0028.0030); In addition, all selected operators are defined in 2D, thus not requiring images adjacent to the one analyzed. nor the data of spacing between cuts; consequently the second and third operators used for the extrusion:! cciñn ele rangs have 3 parameters · d jump of pixels (equal to) And the vertical and horizontal distances between contiguous pixels (values of the mentioned DICOM label).
e) Apply the cl:; sific ~ tlor (4 16). The input of this module is made up of the matrix of the feature vectors generated in the previous stage. This step includes the application of the classifier described in the explanation of the invention for the first embodiment of the invention. The result is in the image being processed. d) Rectify the classification (418) given by the initial mask of the canal, using the following methods: 1) Perform a morphological pattern on the image, consisting of reiterating a dilation followed by erosion. a maximum of 20 iterations using circular structuring elements whose radius depends on the size of the reference region (RR) according to the following criteria: the radius (measured in pixels) is worth rl if size (RR) <a; R2 if a ::: size (RR) <b; r3 if b: 5. size (RR) <c; and r4 if c: 5. size (RR), where rl, r2, r3 and r4 take values close to, 2, 3 and 4 respectively, and a, b and e take values around 300, 600 and 900. respectively. 2) Apply a morphological closure to the image with a maximum of 20 iterations and with structuring elements formed by the 8 3x3 matrices shown in fig. 17. 3) When the set of pixels classified as "conduit ·" forms several unrelated regions, then all the pedals belonging to the smaller regions are reclassified as "background", retaining the classification as "conduit" only for those pixels that belong to the region that contains the largest number of pixels classified as "conduit", and if there is more than one region with the maximum size, the one with its centroid closest to that of the region of interest 4) Reclassify as "background "the pixels that simultaneously meet these three conditions: i) the coordinate of its pixel feature vector generated by the nonlinear operator has the value 0.5. ii) the pixel is on the periphery of the duct; and 3) the coordinate generated by the directional derivative has a positive value. This stage of rectification by knowledge a priori is repeated twice. putting in each one the four steps in different order to minimize the interactions that result from the order in which
different algorithms are applied. e) Control over the number of iterations (420). In this stage, the sum of quadratic differences between the corresponding elements of the two masks was adopted as an index of similarity between two iterations; and as criteria to stop the iterations, that the index is null. If the iteration n: r.: Nien te.ninada is the first one, it is obligatory to make another iteration to obtain a reference with which to compare the segmentation obtained and quantify the similarity or if - being subsequent - the segmf: nulation obtained is not equal to no segmentation obtained before, then another iteration is performed and the method continues with a new iteration, reassessing the specific parameters executed: using step (a) to evaluate and save the values of the ones to start another iteration: from the second iteration, inclusive, if the segmentation obtained is equal to some segmentation obtained before (which may not be the preceding one) then the segmentation obtained as the definitive one for that cut is adopted and the execution continues in the next stage ({) to conclude the segmentation of the current image and continue the segmentation process with the other cuts.
f) Generate the binary channel mask (422). Once the segmentation of the RO1 is finished, the mask of the definitive channel for the analyzed cut is generated, creating a binary matrix that in the pixels! Cs corresponding to the ROl have the values O (background) or 1 (channel) resulting from the segmentation of the ROl, and in the rest it has the value corresponding to the fund
(fig. Ile). Thus, the result of this stage is that the entire image of the cut (and not just the ROl)it is segmented into only two regions that represent the channel and the bottom and that are definedby the definitive mask of the canal (fig. 11 d).
g) Control over the number of segmented images (424): determine if they remainimages pending segmentation; if they remain, return to stage b); if not. continue on stage h)next.h) Generate the channel hypermask (426): the channel mask obtained in each cutit is incorporated into the binary hypermatrix B that rccoge the masks of the channel cn all the cuts;additionally it is shown to the user for tracking and recorded in pemlanentes media andreadable by computer for later use by modeling, measurement andmorphometry
4. MODELING THE NASO LAGRIMAL CHANNEL (206)
4.1. Surface model of the channel (506). The surface model is obtained from the binary hypermatrix 8 that contains the sectional masks of the channel by interpolation with triangles of the contours of the sections of the channel in the successive cuts, based on the three-dimensional metric deduced from the tomographic volume meladatos , softening the surface obtained in the interpolation and finally selecting the isosurface as the channel model whose volume is equal to that obtained in the direct solid model. The numerical structure of faces and vertices thus obtained represents the surface model of the channel: fig. 16 offers two views of the model obtained for the two channels of subject number 17. 4.2. Modeling of the axis of the channel (5/0). 4.2.1. Sequence (the axial Wreo of the channel axis (5/2). The sampling sequence adopted
{CI. and; ..... cp J. is the one given by the centroidt: s of the binary masks of the channel in the successive axial cuts in which it appears.
4.2.2. Axis!! of interpolation of the callal (5/4). Three different interpolation axes are calculated: axial. straight and polygonal.
The axial axis is calculated as the straight segment parallel to the axis of the TAC and between the initial cut to the end of the channel; It is represented by two real numbers that correspond to the z coordinates of the first and last cut medial planes of the channel.
The recl axis () is calculated as the straight segment that joins the centers of the channel in its two extreme cuts and is represented by two tcrnas of real numbers that represent the coordinates of the centroids of the channel in the two extreme cuts of the channel.
The polygonal axis is calculated as the polygonal line that joins the centroids of the axial sections of the channel in all the cuts that contain it.
4.2.3. Eje.'i de aju.'ile del eallal (516). They are calculated by adjusting polynomials of degree 2, 3 and 4 by the method of least squares and I: each polynomial model is represented by the corresponding equation created with the coefficients obtained in the adjustment.
The polynomial model SJ = (x (t) and (t) z (t)) of the axis of a nasolacrimal canal of subject number 26, expressed in parametric form and in the physical coordinate system. is shown in fig. 14a: Figures 14b and 14c illustrate the procedure for adjusting said polynomial model in the sagittal and coronal planes
5. CHANNEL MEASUREMENTS (208)
This stage has the objective of making measurements in an automated way of various physical dimensions of the canal and its sections from its segmentation in the tomography and from the tomography mctadalos. 5.1. 3-D measurements (602)
5.1.1. Volumell of eallal based on seeciolle.'i axiale.'i (604). In this embodiment we do not differentiate primary and secondary voxels, since they coincide when there are no jumps or overlaps between successive cuts. The volume of the voxels has been calculated by multiplying the dimensions of their base from the OlCOM label (0028.0030) by the difference in the third coordinates of "image position (patient)" (0020.0032). The volume of the bone channel for all the subjects of this sample has been calculated by adding the volumes of the voxels that make up the direct solid model of the channel, corrected by the wedge effect at both ends.
5. /.1. Volumell of the eallal base it on the model of! .- uperfh: ie (606) The volume of the channel has also been calculated by numerical methods applied to the surface model, applying the wedge effect correction by the procedure indicated in the section
of explanation of the invention.
5.1.3. Lon¡: illIll llel CUllul (608). In the embodiment described herein, the length is automatically measured by three methods, which name the magnitude obtained. The axial length is obtained by stopping the Euclidean distance between the outer planes of the two extreme cuts of the channel. The length between ends, such as the distance between the centers of the outer planes of the end cuts of the channel. The polygonal length, adding the lengths of the segments that make up the polygonal axis. The polynomial lengths are calculated as the length of the curve corresponding to each polynomial and between the outer planes of the extreme cuts of the channel; The calculation is done by numerical methods with error <0.1 mm. 5.2. 2-D measurements (610)
5.2.1. Section section "to the axial of the channel (612). In each cut this section is measured automatically by counting the number of cut-off axes assigned to the channel by the automatic segmentation program and multiplying it by the area common to all and obtained by multiplying the two There are things given by the OICOM label (0028.0030).
5.2.2. Orthogonal section of the o! Thymus (61-1). It is obtained for each point of the axis sampling sequence. To obtain it, the axis orientation is first calculated at each of these points: for the polygonal axis. the orientation of the axis is calculated by averaging the orientation of the two segments of the polygonal that influence each vertex, considered as free vectors; for the polynomial axis the orientation of the axis at said points is calculated by its tangent at said point and this is calculated from the functions derived from the three components of the parametric expression of the axis (without resorting to numerical approximations); The derived cells are obtained by means of dcrivation rules that are applied automatically during the execution of the program. The orientation of the tangent thus obtained in each case defines the orientation of the plane orthogonal to the axis. Finally, the estimated orthogonal area is obtained by projecting the aforementioned axial area on said plane.
5.2.3. Orthogonal area! (616). It represents the area of the channel in the orthogonal section to the corresponding axis. It is calculated at each point in the sampling sequence of each axis. The result depends on the orientation at which point the axis is used and which determines the orientation of the plane of the cut. the orientation of the axial and inter-axis axes is given by that of the axis itself: that of the polygonal axis at a point p of the sampling sequence is obtained as the vector sum of the segments of the polygonal that converge on p taken as free vectors ; that of each polynomial axis is obtained by calculating the tangent to the curve of the polynomial at point p using its derivative: the orientation of the unit vector in the direction thus obtained is determined by itself the orientation of the new cutting plane, as a vector normal to that plane. The image of the cut in said plane is obtained by applying bicubic interpolation techniques applied to the
images of the original tomography. Finally, the image of the orthogonal cut thus obtained is
apply the automatic segmentation method described above and the mask is obtained
binary channel in the image. The orthogonal sectional area is obtained by multiplying the size of
a pixel given in the metadata by the number of pixels that indicate the channel in the mentioned mask
5 binary of the cond ucto.
The minimum, axial or normal sectional area. represents the minimum value that reaches that area at
along the entire channel. The depth of the minimum area. axial or normal, is the shaft length
measured from the beginning of the channel in the orbit to the point where the corresponding area reaches
Its minimum value.
the
5.2.4. .Axial and orthogonal orientation of the channel section (618). They are calculated by determining
the orientation of the maximum diameter of the section of the channel in the axial section and in the orthogon
respectively. They are expressed as angles in degrees. placing the angular origin in the pane
bottom of the vertical direction of the image; and the positive sense, to the patient's left
fifteen (looking at the cut. the angular origin is down and the angle increases in the opposite direction to
clockwise movement).
5.3. Measurements 1-0 (620)
5.3.1. Minimum and maximum sectional diameter length in axial section (622), which is calculated
twenty in axial images at points in the sampling sequence and is obtained in coordinates
Physics Calling as the size of the pi: <he who gives the DICOM label (0028,0030) and
using standard techniques adapted to obtain measurements between the outer edges of the
pixels instead of between their centers.
5.3.2. Minimum and maximum diameter length in orthogonal cut (624). Dem for
2S minimum distances
6. MORPHOMETRY (210)
This stage (fig. 7) has the objective of obtaining descriptive indexes of the channel's torma and its
sections from the information provided by the numerical models of the channel and its
30 axis obtained in the previous stage.
6.1. Morphometry of the sections of the callal (702) following the procedures indicated in
The description of the invention is cah: ulan:
6.1.1. Momtmlm. · (70-1).
35 6.1.2. Eccentricity in the section of the cimal (706): axial and orthogonal. 6.1.3. Elongation (708).
6.1.4. Circularity (710). 6.2. Conduct shaft morphometry (714)
5 6.2.1. index ... in points of ej e of the channel (716). Using the methods described above, the following indices are calculated. for axial and orthogonal cuts. for the polygonal axis and three polynomial axes. and in all cases in the physical coordinate system established by the metadata: orientation .. curvature and torsion of the channel at all points of its sampling sequence.
10 6.2.2. sect indices: iolJale ...) 'pWlIlIales along the ej e of the channel (7/8). The eccentricity and elongation are calculated for all points of the sampling sequence of each channel, its maximums and minimums are determined, the position of the axis where they are located, the distance they are from the ends of the channel by measuring it at along it, and the overall shape of the canal (double funnel. conical. etc.).
15 For all the specific and sectional variables obtained, their statistics are calculated (mean, maximum, minimum, range and standard deviation) and the graph of each variable has been generated based on the axis of the TAC and based on the axis of the channel
6.2.3. Pllnlm location; singular axis (720). The depth has been determined, starting from the orbit, to which the maximum and minimum of the different areas are located.
20 sectional and curvature and torsion of the channel axis, calculated from the measurements of these variables made t: n all the points of the sampling sequence.
7. SHOW, RECORD AND EXPLOIT THE RESULTS (212)
Once the automated analysis is completed, all results are recorded on all media.
25 obtained and made accessible to the user through graphic screens and printers. For the observation of the models interactive means based on the mouse are used that allow to observe the model from the perspective: and in the lighting conditions chosen by the user. as illustrated in figure 16.
30 The numerical result files with the values obtained for all the described variables were recorded. the format is adapted to make them imported by visualization and statistical analysis programs and transferred to other computers through the network to carry out the visualization and analysis of the results by remote users. A physical copy of the numerical model of the nasocriminal channels is also generated by means of a
35 JO printer.
WAY THAT THE INVENTION IS SUSCEPTIBLE OF INDUSTRIAL APPLICATION
The routes currently planned for industrial application of this invention include: 1) implementation of the method through software programs in medical workstations. for diagnostic and surgical planning assistance, providing accurate representations of the tear ducts; 2) quantitative study of the temporal evolution of the characteristics of the ducts; 3) implementation of a web server capable of receiving a lmomographic study generated in another remote equipment, perform the segmentation and the metric and morphological analysis of the anatomical structures and deliver the models and other results of the study to the applicant of the study through networks local or inter-network telematics or any other means and in the necessary format (3 D images, reports. etc.): 4) incorporation of the methods of the invention to tomognitic devices to guide the process of realization and control of the tomography: 5 ) Endoscopy guidance assistance system; 6) Incorporation into the trajectory planning system and guided exploration or surgical robots; 7) production of tear duct models for teacher use; 8) Production of customized implants from the models generated by the invention; 9) operation of a remote service for the analysis and modeling of lacrimal canals or other related ducts of people and animals through the internet, where the user connects via internet to the site of the company providing the service, transfers the TAC containing the conduit to analyze together with the request to make it. the provider receives it on their server, and it transfers it to the device that performs the three-dimensional study and which has been developed as an embodiment of this invention. It generates the results and transmits them to the applicant through the internet or by any other means.
REFERENCES !. Bookstein FL! 991. Morphometric 10015 for landmark data: geometry and biology Cambridge: Cambridge University Press.
2. Bookstein FL 1997. Landmark m (: thods for forms without landmarks: Localizing group differences in outline shape. Med I mage Anal 1: 225 · 243.
3. Brett M. Hanke M and Larson E. Defining Ihe D1COM orientation hnp: llnipy .orglnibabe I / dicom / dico m _ orientation. hlml; 4.9.201 5.
Four. Clunie, D. Slice spacing of D1COM series. http://public.kitware.com/pipennail/insighlusers/2005-September/0 14711 .html
5. DICOM PS3.3 2015b Information Object Definitions. hnp: llmed ica l. nem a. orgl stan dard. btm l.
6. Estes JL 1. Tsiouris AJ. Christos P.J, Lelli Gl Three-dimensional volumetric assessment
of the naso lacrimal duct in patients with obstruction. Ophthal Plast Reconstr Surg. 2015 May-Jun; 31 (3): 211-4.
7. FelJer, W. An introduction 10 probability theot) 'and its appJications. Vol. JI. 2nd ed. Wiley, 1971
8. Ferson, S., FJ Rohlf, and RK Koehn. 1985. Meas.uring shapc variation of rwodimensional Ollllines. Systematic Zoology. 34: 59-68.
9. Gonzalez R.e. & Woods R E ·· Digitallmage Processing "2nd Ed .. Prentice-f-Iall. New Jersey. 2002
10. Hu M K "Visual Midwife Recognition by Moment Invariants". IRE Transactions on ioformation theory, vollT-8, Page 179.1962
eleven. Jang SW, Seo YJ, Yoo YS, Kim YS.Haahr M., Computed tomographic image analysis based on FEM performance comparison of segmentation 00 knee joint reconstruction. ScieotificWorldJournaL 2014; 2014: 235858. doi: 10.1155 / 2014/235858. Epub 2014 Nov 27
12. Haralick R. M "1974 A Measure fforCircularity of Digital Figures.
13. Huertas A and Medioni G. Detection oflntensity Changes with Subpixel Accuracy Using Laplacian-Gaussian Masks. IEEE Transactions on Pattem Analysis and Machine Intelligence. VOL. PAMI-8. NO. 5. September 1986
14. Klingenberg ep. 1996. Multivariate allomett) '. Ln: Marcus LF .. CoOli M .. Loy A., Naylor GJ., Slice DE. (eds.), Advances in Morphometrics, NATO ASI, Series A: Life Sciences, Vol. 284, pp. 2349. Plenum Press, NY.
fifteen. Klingenberg CP • McJntyre GS. 1998. Geometric Morphometrics of Developmental Instability: Analyzing Panems of Fluctuating Asymmett) 'with Procrusles methods. Evolution 52: 1363-1375.
16. Kroon O J. Segmentarian al / he Mandibular Canal in Cane -Beam CT Data. PHO thesis. An iv. Twente Enschede, 2011).
17. Levine M D Vision In Man And Machine. Chapter 10. "Shape". Page 480. 1985
18. Maatman G. High-Resolulion Computed Tomography of the Paranasal Sinuses and Pharynx and Relaled Regioos Springer Science & Business Media, 1986.
19. Montero and Bribiesca 2009 Intemational MathcmaticaJ Forum, 4, 2009, no. 27,1305 1335 Slatc oflhe Art of Company and Circularity Measurcs
twenty. Pyun J-H, Lim Y-J, Kim M-J. Ahn S-J, Kim J. Position of the mental foramen on panoramic radiographs and ils clause to the Horizontal course of the mandibular canal: a computed tomographic analysis. Clio Orallmpl Res. 24, 2013. 890-895.
twenty-one. Pran William K .. Digitallmage Processing, New York, Joho Wiley & Sons, Ine .. 199 1,
p.634.
22 Ramcy NA. Hoang JK. Richard MJ. CT multidetector of nasolacrimal canal morphology: normal variation by age. gender and raee. Ophthal Plast Reeonstr Surg. 2013 NovOec; 29 (6): 475-80.
2. 3. Rosenfeld A. Compact figures in digital pictures. IEEE Trans.Systems. Man and Cybemeties. 4: 221-223, 1974.
24. Takahashi Y, Nakata K. Miyazaki H, Ichinose A, Kakizaki H. Comparison of bony nasolacrimal canal narrowing wiith or withoul primary acquired nasolacrimal duct obstruelion in a Japanese population. OphthaJ Plast Reconstr Surg. 2014 SepOct; 30 (5): 434-8
25. S. Tangaro.N. Loving, M. Boccardi, S. Bruno, A. Chincarini, G. Ferraro, G.B. Frisoni, R. Maglictta, A. Redolfi, L. Rei, A. Taleo, R. Bellotti Automated voxel-by-voxel tissue classification for hippoeampal segrnentation: Mcthods and vaJidation. Physica Medica: European Journal of Medical Physics. Volume 30. Issue 8,878-887 2014
26. Yong AM lhao OB, Siew Se. Goh PS, Liao J. Amrith S. Assessment of bony nasolaerimal parameters among Asians. Ophthal Plast Reconstr Surg. 2014 JulAug: 30 (4): 322-7.
27.2unic. J .. Hirota K. Measuring Shalpe Cireularity Progress in Pattem Recognition, Image Analysis and Applications. Lecrurc: Notes in Computer Science Yolume 5197, 2008, pp 94-101

权利要求:
Claims (1)
[1]
l. Method to quantify. departing dI! a CT scan the metric and morphological characteristics of a duct of irregular fonna and of heterogeneous walls and content, and that is characterized by automating the pro: duct segmentation process by applying to each tomography image an iterative processing that in each iteration includes the four stages following:
(to) calculate the trait vector associated with each pixel of the region of interest (ROl). vector that will have at least dimension 2, each of whose coordinates is the value resulting from applying in said pixel the corresponding feature extraction operator, calculated from the value of the pixels and of the general parameters and the specific parameters of each iteration which are recalculated at the beginning of it:
(b) classify each pixel of the RDI by labeling it with the value I or O • representative respectively of "duct" and "background", resulting from applying a binary classifier to the feature vector of said pixel: and creating the matrix of the initial duct mask with the values thus obtained in the pixels of the ROl and with Oen the others:
(C) rectify the classification given by the initial duct mask by inverting the binary value of its elements corresponding to the RDI when they do not conform to the pre-existing knowledge about the nature of the substances constituting the tomographed material and on the shape of the contact surface between the duct and the bottom; The matrix thus obtained is the rectified duct mask:
(d) iteration control: compare the current rectified duct mask with the rectified masks obtained in all previous iterations in the same image; if the pre-established similarity criterion is not satisfied in any comparison. then continue on the stage
(to) to make an additional iteration: otherwise add the current rectified mask. as a definitive duct mask in the analyzed image, to those obtained previously. to build with them all the binary hypernatrix that defines the segmentation of the duct.
The molasses according to claim 1 also includes the previous Napa to delimit the region of interest (ROl) for each and every one of the images. region that will include the entire section of the duct and an additional part of the duct: the surrounding area. being able to understand the whole image; said region will be defined by a mask that will have the value l in the position of the pixels that pc: belongs (; in the region of interest and Ocn the rest; the delimitation of the ROl can be done by automatic means. manual or assisted by computer.
The method according to claim 2 in td that the delimitation of the computer-assisted RDI includes one or more of the following aids: zoom with inerpolation, smoothing of the image, superimposing reference lines on the image. and visualize the image with a level in a range of background values close to the gray levels inside the duct and with a window around a third of the range of background values.
4 The method according to claim 3 in which the delimitation of the computer-assisted ROl includes zooming around x8 with interpolation of the nearest neighbor, smoothing the image by a low-pass Gaussian filtering, superimposing reference lines on the image that They include crossing the dellaplaciano by zero or an iso-contour close to 175, and visualize the image with a level close to 300 and a window around 1,000.
The method according to claim 2 wherein the delimitation of the ROl is performed automatically.
The method according to claim 1 wherein the general parameters comprise: a) the numbers of the first and last image of the tomographic series that intersect the ohjeto; b) the nature of the similarity index between two successive segmentations; c) the cut-off value required for said index to stop the segmentation iterations with each image; d) the criteria for determining at time of execution the shape and size quc should have the structuring elements in the operations of mathematical morphology; and e) the ranges of the values provided in the tomography for each type of substance that integrates the tomographed object, derived from the preexisting knowledge.
The method according to claim 1 wherein the specific parameters for the current iteration comprise: a) determining the region of reference, which for the first iteration of segmentation in a tomographic image is the ROl. and for the second and subsequent ones it is determined by the rectified duct mask obtained in the preceding iteration; b) determine the area of the reference region in the physical coordinate system. using for this the information about the pixel size contained in the metadata; e) calculate the physical coordinates (c, CI) of the central reference point, given by the cenlroid of the reference region; and d) determine the shape and size of the structuring elements to be applied in the operations of mntematic morphology - erosion and dilation - calculated according to the size of the reference region and the criteria established in the general parameters.
The method according to claim 1, whose step of calculating the feature vector also includes as a stage prior to feature extraction reducing image noise by means of a median filter. or a Gaussian low-pass filtering or the convolution of the image with an optimized matrix to eliminate the type of electronic noise from the tomography.
The method according to any of the preceding claims in which the feature vector includes as one of its coordinates the value calculated in each pixel by a non-linear operator. 2-D or 3D. whose parameters are dynamically recalculated for each iteration of segmentation; and as another of its coordinates the value adlJptado in said pixel by a discrete differential operator calculated using the physical coordinate system determined by the mt: tadata.
The method according to claim 9 wherein the non-linear feature extraction operator consists of a threshold with at least two thresholds. I1 and (J. whose values depend simultaneously on the general parameters set a priori and on a central tendency statistic of the image values in a local environment of the central reference point; said operator calculates for each pixel the value of the trait in function of the value of said pixel. X,}. assigning the value "1" when X,} <(/, the value "0" when Xi)> tl Y the value "0.5" otherwise.
11 The method according to claim 10 wherein the umh rales ti and 12 are calculated as follows: I1 = min (Cma '. Max (Cmm. Mu)] and t; & (Bmm + (1) / 2, where Mu is a central trend statistic of the values of a submatrix nxn of the image centered on the central reference point, taking n a value around 3; B "" n represents the minimum anticipated value for pixels not belonging to the conduit .and C nun and Cm. «represent respectively the lower and upper dimensions for the range of the maximum value expected in each cut for the pixels belonging to the" duct "; the values of CMtu, c '", n and 8mm will have been established by an expert depending on the nature of the tomography duct and the type of tomography used.
The method according to claim 11 wherein Cmn, CIII <U and 8 "," 1 take values in the intervals whose approximate limits are respectively 30-80.80-120 and 300-700.
The method according to any one of the preceding claims in which one of the feature extraction operators consists of an active contour that a) is initialized with the same form 'j center as the reference region' j with a size equal to one fraction of it; b) the contour evolves towards the periphery of the channel attracted by a function of the partial value derivatives
of the pixels until it reaches its f! stationary state; and c) the nonlinear operator assigns the value I or O to the trait according to whether the pixel is located inside or outside - respectively - of the stationary contour, and 0.5 otherwise.
The method according to any of the preceding claims in which one of the feature extraction operators consists of an operator based on level or basin sets ("watershed") that the) applies an algorithm to the RlDI of segmentation using the technique of basin-level sets ("watershed"); and 2) the operator assigns the trait to the value 1. It is 0.5 that the pixel is located in the greater of the interior regions, in the greater of the peripheral regions, or in another, respectively.
15 The method according to any of the preceding claims in which the differential feature extraction operator consists of a Laplacian operator who in each pixel of the ROl assigns the discrete Laplacian value of the image values on the said pixel to the corresponding feature.
The method according to any one of the preceding claims wherein one of the feature extraction operators consists of the directional derivative appeller calculated by means of the discrete derivative in the direction of an identifiable reference point located inside the conduit.
17 The method according to claim J 6 wherein the identifiable reference point is the centroid of the reference region in each iteration.
The method according to any of the preceding claims whose classification step includes a deterministic decision tree that labels each pixel of the ROl as a "conduit" or as a "background" based on the values of its feature vector and comprising the following steps i) label as "conduit" the pixels in which the feature vector takes the value 1 in the coordinate generated by the nonlinear operator; label as "background" those in which said coordinate is worth O; and label those in which 0.5 is worth as "conduit" or as '· background' according to whether the coordinate from the differential operator is positive or negative, respectively; ii) classify the rest of pixels with 0.5 in its first coordinate as "conduit" or "background", when in a small one or in its own there is a majority of pixels already classified as "conduit" or as "background", respectively; and when both ~ proportions ~ are equal then it will be classified as "background" when the average value of [pixels in the
cited environment is closer to the higher threshold and as a "conduit" in another case.
The method according to any of the preceding claims in which the rectification step comprises one or more of the following operations: 1) applying a morphological closure, consisting of reiterating, with a limit on the number of iterations, a dilation followed by erosion . with circular structuring elements whose radius is in inverse relation to the maximum curvature admitted in the separation surface of the "figure" and the "bottom" and in direct relation to the size of the reference region; 2) apply a morphological closure with a limit on the number of iterations and with fibrillar and small size strruciant elements; 3) when the set of pixels classified as "I: wavy" forms several unconnected regions then all pixels belonging to the smaller regions are reclassified as "background", retaining the classification as "conduit" only for those pixels belonging to the region that contains the largest number of pixels classified as "conduit", and if there were more than one region with the maximum size. the one with its centroid closest to that of the region of interest; 4) reclassify as · · background "105 pixels that simultaneously meet the following conditions: i) the coordinate of its feature vector generated by the nonlinear operator has the value 0.5, ii) is n at the periphery of the condlJcto.
The method according to any one of claims I to 18, wherein the rectification step comprises 1) Performing a morphological closure in the image, consisting of reiterating a dilation followed by erosion. a maximum of 20 iterations using circular structuring elements whose radius depends on the size of the reference region (RR) according to the following criteria: the radius (measured in pixels) is given if size (RR) <a; R2 if a.5: size (RR) <b; r3 if b.5: size (RR) <c; and r4 if c.::5. size (RR), where r1. r2, r3 and r4 take values close to 1,2.3 and 4 respectively, and a. b and c take values around 300, 600 and 900, respectively. 2) Apply a morphological closure to the image with a maximum of 20 iterations and with structuring elements formed by the 8 3x3 matrices shown in fig. 17.3) When the set of pixels classified as "conduit" contains several unrelated regions, then the pixels belonging to the smaller regions are reclassified as "background", retaining the classification as "conduit" only for those pixels is that they belong to the region that contains the largest number of pixels classified as "conduit". and if there were more than one region with the maximum size, the one with its centroid closest to that of the region of interest. 4) Reclassify as "background" the pixels that simultaneously meet these three conditions: i) the coordinate of its pixel feature vector generated by the nonlinear operator has the value 0.5. ii) the pixel is at pt: duct rifery; and 3) the coordinate generated by the directional derivative has a positive value.
21 The method according to any of the preceding claims which further includes the method for obtaining the image of the axial cut at an arbitrary point Zo of the z axis of the tomography. along with their metadata. and the corresponding binary mask of the duct section in said image. and that includes the following steps: a) if the value of Zo is equal to the z coordinate of any of the cuts of the original tomography, then the desired result will be formed by the image. the metadata and binary mask corresponding to said cut; b) in another case, calculate the multiplanar reconstruction while maintaining the orientation and spacing of pixels and cuts, but matching Z {J the z coordinate of one of the cuts; apply to the multiplanar reconstruction the preprocessing and segmentation method; the result sought is given by the image of the cut in Zn resulting from the multiplanar reconstruction, along with its metadata. and the corresponding binary mask in the binary hypernatrix resulting from segmentation
22 The method according to any of the preceding claims which further includes the method for calculating the direct solid model. which is the geometric object formed by the set of "secondary voxels" corresponding to IIOS pixels of all the cuts that have been classified as "conduit" in the segmentation process, meaning "secondary voxel" corresponding to a pixel the solid geometric body defined in the physical space whose centroid and size of the base in 3D physical coordinates coincide with those of the pixel and whose height is given by the spacing between the centers of successive cuts.
The method according to claims 1. or 22 which further includes the method for calculating the volume of the duct in the physical coordinate system that is obtained as the sum of the volumes of all the "secondary voxels" that form the direct solid model.
The method according to claims 1 or 22 which further includes methods for generating: a) the surface model that is obtained from the binary hypermatrix B containing the sectional masks of the duct interpolating the contours of the duct sections by means of a surface in the successive cuts, based on the three-dimensional metric deduced from the metadata of the tomographic volume. softening the surface resulting from the interpolation, and selecting as a conduit model the isosurface that defines a volume equivalent to that based on the direct solid model. The numerical structure of faces and vertices thus obtained represents the surface model of the conduit; or b) the solid fill model, which is obtained as the closed set of the space delimited points for each surface model.
The method according to claim 24, which further includes the methods for calculating the volume of the duct delimited by the surface model or the solid filling model using the numerical or finite element integration techniques.
26 The method according to claims 23 or 25 which further comprises: a) calculating the error by wedge effect that by excess or by default h: a has been included in the value of the volume obtained and (b) apply the corresponding correction to the volume obtained.
27 The method according to claims 23 or 25 which further comprises, when the conduit is open at one or both ends, apply a correction by wedge effect to the volume obtained which for each end is calculated thus · a) determine the outer plane of the extreme axial cone! '> that contain the duct; b) determine the plane most distal and orthogonal to the axis of the duct whose intersection with the surface that delimits the duct is a closed curve; e) calculate the volume of the wedge between the outer flat cylindrical. the aforementioned more distal plane and the outer surface of the duct; d) subtract said correction to the volume of the rest of the duct if said wedge was already contained in it or add it otherwise.
28 The method according to any of claims 1,9,10,11,12,18,19,20,22 OR 24 which further includes the method for gt = denying a plurality of sampling sequences including a) axial sampling sequence , obtained as a succession of the centroids, expressed in the physical coordinate system, of the axial sections of the duct in successive axial cuts; b) Sequence of skeletal death, obtained as a succession of the intersection points, expressed in the physical coordinate system, of the medial planes of the axial slices of the tomography with the symmetry axes, in the axial direction of CT . of the voxcles belonging to the resulting pruned skeleton: from: applying to the direct solid model any of the methods of thinning or skeletonizing images; or c) the derived sampling sequences obtained as moving weighted moving average or other smoothing methods applied to any other sampling sequence.
29 The method according to claim 28 which further includes methods for generating a multiplicity of models of the axis of the conducum that include: a) the interpolation axes. which are calculated by interpolating the points of the primary sampling sequence, OR of the derived sequence. by any interpolation method; or b) the adjustment axes, which are calculated by adjusting functions
vector with values in IRJ to the points of the sampling sequence, or of the derived sequence, by least squares or any other method of adjustment;
The method according to any one of claims 1,9, 10, 1, 12, 18, 19, 20 or 22 which further includes the method for generating the skeleton axis by applying any method of thinning to the direct solid model.
The method according to claim 28 also includes calculating a) the polygonal axis obtained in the following steps: (i) determining the polygonal joining the points of the sampling sequence of the axis and (ii) prolonging each of its extreme segments in an amount equal to half its length, this polygonal being enlarged the polygonal axis of the duct; or b) the right axis, which is given by the rectilinear segment that joins the two ends of the polygonal axis; or c) the interpolation axis sync. or Bezier, or NU'RBS (non-uniform rational B-splines), or spline given by the parametric curve obtained by applying the sine or Bezicr or NURBS or splline interpolation method to the sampling sequence, respectively.
32 The method according to claim 28 which further includes calculating the axis of the conduit given by a vector function whose components are polynomials of the smallest degree with which an acceptable adjustment to the points of the axis sampling sequence is obtained.
33 The method st: gun any dt: the n: claims: from 29 to 32 which also includes the methods of differential and numerical calculation for carrying out the analysis of the axis of the duct and comprising evaluating: a) the length of the conduit, which is obtained by calculating the length of the curve that represents the axis adopted as a model and between the planes orthogonal to the conduit that delimits it at its two ends; b) the orientation at each point of the axis and the average orientation at the points of the sampling sequence; c) instantaneous, average and maximum curvatures; d) instant torsions. average and maximum; e) the location of the points of maximum and minimum curvature and torsion; f) the Euclidean distances from said singular points to each end of the axis; and g) the distances along the axis between said singular points and each end of the axis.
34 The method according to any of claims 29 to 32 which further includes the method for obtaining the image of the orthogonal cut to the axis of the duct at an arbitrary point p thereof, together with the corresponding metadata, and the corresponding orthogonal binary mask, and that
It includes the following stages: a) calculate the tangent in p to the mathematical function adopted as
shaft model; b) determine the equation of the plane perpendicular to the axis of the duct in p; c) calculate the multiplanar reconstruction containing dic ho plane perpendicular to the p-axis, together with the corresponding mctadata, and obtain the orthogonal sectional image of the duct at point p along with its metadata. what will it be equal to the image of the cut corresponding to said point in said multiplanar reconstruction. together with its metadata: d) apply to the reconstruction multiplanar the segmentation method of this invention and obtain the orthogonal binary duct mask in the cut perpendicular to its axis at point p that will be equal to the binary matrix corresponding to point p included in the binary hypernatrix obtained.
The method according to claims: 21 or 34, which also includes the methods for obtaining a plurality of measurements 20 of the section of the pipeline by an arbitrary point p of the axis thereof, each of which will be axial or orthclgonal according to the section on which it is obtained is axial or orthogonal respectively. measurements calculated all of them in the image coordinate system or in the physical coordinate system determined by the metadata. and that include one or more of the following measurements: a) sectional area of the duct in the cut by point p, obtained by adding the areas of the pixels marked "duct" in the corresponding binary duct mask, the area of a pixel equal to the product of its base by its height expressed in physical coordinates by tomography meladatos: b) estimated area of the duct section in the orthogonal plane, obtained by multiplying the corresponding area in the axial section by the cosine of the angle between the axis of the lomography and that of the duct; e) maximum or minimum sectional diameter length, calculated as the maximum or minimum distance - respectively - between the outer edges of all the pairs that can form the pixels located on the periphery of [a section of the duct; d) orientation of the duct section, given by the orientation of its maximum diameter: or e) moments of the section. standards. centrals, unnalized and invariant centrals. calculated by applying the universally accepted definitions; f) eccentricity of the section. calculated as the ratio between the distance between the foci and the length of the major axis of an ellipse with nominal second central moments equivalent to those of the section of the duct in the image of the cut, or by any other available procedure; g) elongation of the section of the duct, calculated by subtracting from the unit the quotient between the minor and major sides of the envelope rectangle of the section in the chosen coordinate system. or by any other available method; h) circularity of the s (: ction. calculated as the inverse of 2p multiplied by the quotient between the square of the area and the sum of the variances in the two axes; i) contour line of the section, calculated as the set of points of contact with the bottom of the
pixels that belong to the duct; j) the plurality of indexes of the section that offers the
Calculated geometric morphometry direc [amenle on a contour line of the duct or starting from marks or "Iandmarks" derived from it.
The method according to claims 3 3 or 35 which also includes obtaining statistics of each variable measured according to said claims by: a) repeated carrying out the measurements of said variable at all points P, of the sampling sequence of the axis . or other set of axis points that are of interest: b) the calculation of the statistics of the measurements obtained. which include one or more of the following indices: i) the indices of central tendency and variability, ii) their extreme, minimum and maximum values; iii) determine the points of the axis where these external values are located; and iv) obtain the curves of each variable that express the measured value in ordinates and in abscissa the axis of the tomography or the rectified curve of the axis of the duct.
37 The method according to any of the preceding claims, which also includes the methods for a) calculating the corresponding contour line from each sectional binary or orthogonal binary mask. ~ Axial or orthogonal · respectively · in the coordinate system chosen one. line that as such has a null thickness and that is integrated by the set of contact points between the pixels of the duct and those of the bottom expressed in the physical coordinate system; b) generate from it the marks ("Iandmarks") and the matrices of reference points in physical coordinates established in the metadata of the binary image contour of the conduit in the orthogonal cut to its axis, and rt: prescfltado CI1 la corresponding binary sectional image: and c) calculate a plurality of geometric morphometry indices 20 based on said marks.
38 The method according to any of claims from] to 36 which further includes the methods for: a) obtaining the marks ("Iandmarks") and matrices of reference points, in the physical coordinate system established in the metadata, from from solid 3 · 0 models or from the surface model or from the ~ ie model of the duct; and b) calculate a plurality of 3 O geometric morphometry indices based on said reference marks and matrices.
31; 1 The method according to any of the preceding claims, repeatedly applied to the successive tomographs of the same conduit obtained over time, which also includes the methods for a) calculating evolutionary indices of each variable and b) generating the curves and other formal models representative of the eve, lueion dc said variables over time.
The method according to claim 34 which further includes the creation of [a rectified tomography, following these steps; a) obtain the image of the orthogonal section a [duct axis at all points in the axis sampling sequence; c) position each image thus obtained in a plane perpendicular to the axial axis that [or cut by the z coordinate corresponding to the image and moving it in the indicated plane so that the centroid of the duct section coined with the axis, obtained thus the rectified tomography of the duct; d) record it in permanent memory; e) represent it graphically for visual and morphometric analysis.
The method according to any of the preceding claims, which further comprises methods for local or remote mode a) pennanently register the models obtained;
or b) make them accessible to the user, visualizing the models in virtual reality devices, screens or paper, in their entirety or by cuts; isolated or integrated in a volumetric rendering 11 another representation of its environment; In video mode. static figures graphics. or interactive figures allowing the user to select the observation point. the colors and nature and position of the light sources; or c) print in 2D or 3D of the results chosen by the user; or d) generate digital copies on fixed or removable storage media.
42 The method according to any of the preceding claims, further comprising methods for locating or remotely a) recording in a maneral way the graphic and numerical results obtained; or b) make them accessible to the user. visualizing them on screens, paper or other volatile or permanent physical supports. in its entirety or by palms; in video mode, static figures, graphics. or interactive figures; or e) print the results chosen by the user; or d) generate digital copies on fixed or removable storage media.
The method according to any of the preceding claims, which also includes the methods for modifying the format of the models obtained so that they are processable by 3D printers, numerical control machines or other means of machining or producing objects from numerical models. and capable of generating from the surface model a custom implant made in the type of material indicated by an expert or a physical representation of the duct model in a plurality of scales and materials.
44 The method according to any of the preceding claims which further includes the methods of calculating trajectories for guiding robotic devices, including surgical ones, in cleaning, draining, filling, repairing or surgery of the ducts represented in the
surface model.
45 A programmable storage device that is readable by a computer, not volatile. which tangibly incorporates a program of executable instructions in a computer, a digital signal processor, an application-specific integrated circuit, a microprocessor, a microcontroller, or any other form of programmable hardware, with working memory, media entry. output, internal data storage, and with communications means to import data and export results and able to perform the steps of the method to quantify, from a computerized tomography, the metric and morphological characteristics of an irregularly shaped duct and walls and heterogeneous content, said method being characterized by automating the duct segmentation procedure by applying to each tomography image an iterative processing that in each iteration includes the following four stages:
(to) calculate the associated feature vector; to each pixel of the region of interest (ROl). vector that will have at least dimension 2, each of whose coordinates is the value resulting from applying in said pixel the corresponding feature extraction operator, calculated from the value of the pixels and the general parameters and the specific parameters of each iteration that are recalculated at the beginning of it:
(b) classify each pixel of the ct ROl: labeling it with the value I or O • representative respectively of "duct" and "background", resulting from applying a binary classifier to the feature vector of said pixel; and create the matrix of the initial duct mask with the values thus obtained in the RDI pixels and with O in the others;
(and) rectify the classification given by the initial duct mask by inverting the binary value of its elements corresponding to RO l when they do not conform to pre-existing knowledge about the nature of the substances constituting the tomographed material and about the shape of the contact surface between the duct and the bottom; the matrix thus obtained is the rectified duct mask;
(d) iteration control: compare the current rectified duct mask with the rectified masks obtained in all previous iterations in the same image; If the pre-established similarity criterion is not satisfied in any comparison, then continue on the stage
(to) to do an additional iteration; otherwise, add the current rectified mask, as a definitive duct mask in the analyzed image. to those obtained previously. to build with all eHas the binary hypermatrix that defines the segmentation of the duct.
The device according to claim 45, which also includes the means for carrying out the previous stage of delimiting the region of interest (ROl) for each and every one of the images, a region that will include the entire section of the condom and an additional part of the substance that surrounds it, being able to understand the whole picture; said region will be defined by a mask that will have value 1 at the position of the pixels belonging to the region of interest and O in the rest; The limitation of the ROl can be done by automatic, manual or computer-aided means.
47 The device according to claim 46 wherein the delimitation of the computer-assisted ROl includes one or more of the following aids: interpolation zoom. image smoothing superimpose reference lines on the image, and display the image with a level in a range of background values close to. gray levels inside the duct and with a window around a third of the background value range.
The device according to claim 47, wherein the delimitation of the computer-assisted ROl includes zooming around x8 with interpolation of the nearest neighbor, smoothing the image by a low-pass Gaussian filtering, superimposing reference lines on the image that Inc They include the zero crossing of the Laplacian or an isocontom close to [75, and visualize the image with a level close to 300 and a window that took 1,000.
49 The device according to claim. 46 in which the means for realizing the delimitation of the ROl by an automated procedure are included.
The device according to claim 45 wherein the general parameters comprise: a) the numbers of the first and the last image of the tomographic series that intersect the object; b) the nature of the similarity index between two successive segmentations. c) the cut-off value required at said index to stop the segmentation iterations with each image; d) the criteria to determine at runtime the form and size that the structuring elements must have in the operations of mathematical morphology; and e) [the ranges of values provided in the tomography for each type of substance that integrates the scanned object, derived from pre-existing knowledge.
The device according to claim 45 wherein the specific parameters for the current iteration comprise: a) determining the reference region, which for the first iteration of segmentation in a tomographic image;: I. it is the RD [, and for [a second and following it is determined by the rectified duct mask obtained in the preceding iteration: b) determine the area of the reference region in the physical coordinate system, using for
this is the information about the size of the pixels contained in the metadata; e) calculate the physical coordinates (e, e) of the reference center point I, given by the centroid of the reference region; and d) determine the shape and size of the structuring elements to be applied in the operations of mathematical morphology - erosion and dilatation - calculated according to the size of
5 the reference region and the criteria established in the general parameters.
52 The device according to claim 45. wherein the step of calculating the feature vector further includes as a stage prior to feature extraction reducing image noise by means of a median filter, or a Gaussian low pass filtering or the convocation of the image
10 with an optimized matrix to eliminate the specific noise pipe from the tomography.
The device according to any of claims 45 to 52, wherein the feature vector includes as one of its coordinates the value calculated in each pixel by a non-linear operator. 2-D or 3D, whose parameters are dynamically recalculated for each iteration of
15 segmentation; and as another of its coordinates, the value adopted in said pixel by a discrete differential operator calculated using the physical coordinate system determined by the metadata.
54 The device according to claim 53 wherein the non-linear operator of extraction of
20 features consist of a threshold with at least two thresholds, ti and J2, whose values depend simultaneously on the general parameters set a priori and a central tendency statistic of the image values in a local environment of the central reference point; said operator calculates the value of the feature for each pixel based on the value of said pixel. x '/' assigning the value "1" when x I) <ti, the value "O" when x I}> 11 And the value "0.5" otherwise.
The device according to claim 54 wherein the thresholds t / and 1 / are calculated as follows: tI = min rCmal ". Max (C .. ,,,, M ,,)} Y 11 = (Bm, ·" + ti) / 2, where M¡J is a central tendency statistic of the values of a sub-matrix nxn of the image centered on the central reference point, taking n a value at spine a3; BIft ... represents the minimum value anticipated a priori for
30 pixels not belonging to the conduit, and Cm ", and CntlU respectively represent the lower and upper dimensions for the range of the maximum value expected in each cut for the pixels belonging to the ·· conduit · '; the values of C"' aI " 'Cm ... and B "" "will have been established by an expert depending on the nature of the tomography duct and the type of tomography used.
The device according to claim 55 wherein Cmm, Cmar and BII / m take values in the intervals whose approximate limits are respectively 30-80, 80-120 Y300-700.
The device according to any one: of claims 45 to 56 wherein one of the feature extraction operators consists of an active contome that a) is initialized with the same shape and center as the reference region and with a size equal to a fraction of it: b) the contour evolves towards the periphery of the channel attracted by a function of the partial derivatives of the pixel value until it reaches its steady state; and c) the nonlinear operator assigns the value I or O to the trait according to whether the pixel is located inside or outside respectively of the stationary contour, and 0.5 otherwise.
58 The device according to any of the claims: s from: 45 to 57 in which one of the feature extraction operators consists of an operator based on basin level assemblies ("watershed") that O) applies to the RDI a segmentation algorithm using the technique of level or basin sets ("watershed"); and 2 °) the operator assigns the value 1. 1. 0 0.5 according to the pixel being located in the largest of the interior regions. cn the largest of the peripherals, or in another, respectively.
The device according to any one of claims 45 to 58, wherein the differential feature extraction operator consists of a Laplacian operator who in each pixel of the ROl assigns the corresponding Laplaeian value to the corresponding feature of the values of [image]. on the said pixel.
The device according to any one of claims 45 to 59, wherein one of the feature extraction operators consists of the directional derivative operator calculated by the discrete derivative in the direction of an identifiable reference point located in the interior of the conduit.
61 The device according to claim 60 wherein the identifiable reference point is the center of the reference area in each iteration.
62 The device according to any of claims 45 to 61 whose classification stage includes a deterministic decision tree that labels each pixel of the ROl as a "conduit" or as a "background" based on the values of its feature vector and comprising the following steps i) label as "conduit" the pixels in which the feature vector takes the value I in the coordinate generated by the non-linear operator; label as "background" those in which said coordinate is worth O; and label those in which 0.5 is worth as "conduit" or as "background" according to whether the coordinate from the differential operator is positive or negative, respectively; ji) classify the rest of pixels with 0.5 in its first coordinate as "conduit" or "background", when in a small environment of its own there is a majority of pixels already classified as "conduit" or as "background", respectively: and when both proportions are equal then it will be classified as "background" when the average value of the pixels in said environment is closer to the higher threshold and as "conduit" in another case.
63 The device according to any of claims 45 to 62, wherein the rectification step comprises the means for performing one or more of the following opcrations: 1) applying a morphological closure, consisting of reiterating, with a limit on the number of iterations , a dilation followed by erosion. with the circular structuring elements whose radius is in inverse relation to the maximum curvature admitted in the separation surface of the "figure" and the "bottom" and in direct relation to the size of the reference region; 2) apply a morphological closure with a limit in the number of iterations and with fibrillar and small size structuring elements; 3) when the set of pixels classified as "conduit" forms several unrelated regions then all pixels belonging to the smaller regions are reclassified as "background", retaining the classification as "conduit" only for those pixels belonging to the region quc contains the largest number of pixels classified as "conduit", and if there is more than one region with the maximum size, the one with its centroid closest to that of the region of interest; 4) reclassify as "'background" the pixels that simultaneously meet the following conditions: i) the coordinate of its feature vector generated by the nonlinear operator has the value 0.5, ii) they are on the periphery of the conduit.
64 The device according to any one of claims 45 to 62, wherein the rectification stage comprises the means for: 1) Performing a morphological closure in the image, consisting of removing a continuous expansion of an erosion, a maximum of 20 iterations using circular structuring elements whose radius depends on the size of the reference region (RR) according to the following criteria: the radius (measured in pixels) is worth rI if size (RR) <a: r2 if a'S size (RR) < b: r3 if b.:s: size (RR) <c; and r4 if c.:s: size (RRl. where r 1, r2, r3 and r4 take values close to, 2.3 and 4 respectively. and a. bye take values in volume at 300. 600 and 900. respectively. 2) Apply a morphological closure to the image with a maximum of 20 iterations and with structuring elements formed by the 8 3x3 matrices
collected in fig. 17.3) When the set of pixels classified as "conduit" forms
several unconnected regions are then reclassified as "background" all pixels belonging to the smallest regions, retaining the classification as "conduit" only for those pixels belonging to the region containing the largest number of pixels classified as "conducle", and if there were more than one region with the maximum lamano, the one with its cenlroid closest to that of the region of interest. 4) Reclassify the pixels that meet as background
simultaneously these three conditions: i) the coordinate of its pixel feature vector generated by the nonlinear operator has the value 0.5. ii) the pixel is on the periphery of the duct: and 3) the coordinate generated by the directional derivative has a positive value.
The device according to any of claims 45 to 64 which further includes the means for obtaining the image of the axial cut at an arbitrary point z () of the z axis of the tomography. along with its metadata, and the corresponding binary mask of the duct section in said image. and that includes the following steps: a) if the value of z () is equal to the z coordinate of any of the cuts of the original tomography. Then the desired result is formed by the image. the metadata and binary mask corresponding to said cut; b) in another case, calculate the multiplanar reconstruction while maintaining the orientation and spacing of pixels and cuts, but matching Zn to the z coordinate of one of the cuts; apply to the multiplanar reconstruction the preprocessing and segmentation method; the result sought is given by the image of the z-cut (resulting from the multiplanar reconstruction, along with its metadata, and the corresponding binary mask in the binary hypermalriz resulting from the segmentation.
66 The device according to any of claims 45 to 65 which further includes the means for calculating the direct solid model, which is the geometric object formed by the set of "secondary voxels" corresponding to the pixels of all eons that have been classified as a "conduit" in the segmentation process. "secondary voxel" corresponding to a pixel is the solid geometric body defined in the physical space whose centroid and size of the base in physical coordinates 3 D coincide with those of the pixel and whose height is given by the spacing between the centers of successive cuts .
67 The device according to claims 45 or 66 which further includes the means for calculating the volume of the duct in the physical coordinate system that is obtained as the sum of the vo lumencs of all the '' secondary voxels '' that form the direct solid model .
68 The device according to claims 45 or 66 which also includes the means for generating:
a) the surface model that is obtained from the binary hypermatrix B that contains the sectional masks of the duct interpolating the contours of the duct sections in the successive cuts by means of a surface. based on the three-dimensional metric deduced from the lomographic volume metadata, softening the surface resulting from the interpolation, and selecting the isosurface that defines a volume equivalent to that based on the direct solid model as the conduit model. The numerical structure of faces and vertices thus obtained represents the surface model of the duct; or b) the solid fill model, which is obtained as the closed set of the points of the space delimited by each surface model.
69 The device according to claim 68 which further includes the means for calculating the volume of the duct delimited by the surface model or the solid filling model using the techniques of numerical integration or finite element.
The device according to claims 67 or 69, which further comprises the means for: a) calculating the error due to wedge effect which, by excess or by default, has been included in the value of the volume obtained and (b) apply the corresponding correction to the volume obtained .
The device according to claims 67 or 69, further comprising, when the conduit is open at one or both ends, the means for applying a wedge effect correction to the volume obtained which is calculated for each end as follows: a) detem, inar the outer plane of the extreme axial cuts that contain the duct; b) determine the plane most distal and orthogonal to the axis of the duct whose intersection with the surface that delimits the duct is a closed curve: c) calculate the wedge volume between the said outer plane, the aforementioned more distal plane and the supe : outer surface of the duct; d) subtract said correction to the volume of the rest of the duct if said wedge was already contained in it or add it otherwise.
72 The device according to any of claims 45, 53, 54, 55, 56, 62, 63, 64, 66 or 68 which further includes the means for generating a plurality of sampling sequences including a). of axial axial life. obtained as a succession of the centroids. expressed in the physical coordinate system. of the axial sections of the duct in the successive axial cuts; or b) Sequence of skeletal sampling, obtained as a succession of the intersection points, expressed in the physical coordinate system, of the medial planes of the cuts
axes of the tomography with the axes of their metrics. in the axial direction of TAC, of the voxels
Penetrating the pruned skeleton resulting from applying to the direct solid model any of the methods of thinning or skeletonizing images; or e) the derived sampling sequences obtained as moving average, weighted moving average or other smoothing methods applied to any other sampling sequence.
The device according to claim 72, which further includes the means for generating a multiplicity of models of the axis of the duct including: a) the interpolation axes. which are calculated by interpolating the points of the primary sampling sequence, or the derived sequence. by any interpolation method; or b) the adjustment axes, which are calculated by adjusting vector functions with values in ~: ¡to the points of the sampling sequence, or of the derived sequence, by least squares or any other adjustment method;
74 The device according to any of claims 45,53,54,55,56.62,63,64 or 66 which further includes the means for generating the skeleton axis by applying any method of thinning to the direct solid model.
The device according to claim 72 which further includes the means for calculating a) the polygonal axis obtained in the following steps: (i) determining the polygonal joining the points of the sampling sequence of the axis and (ii) prolonging each one of its end segments in an amount equal to half its length, this polygonal being enlarged the polygonal axis of the duct; or b) the straight axis, which is given by the rectilinear segment that joins the two ends of the polygon axis or e) the interpolation axis sine, or Bezier, or NURBS (non-uniform rational B-splines). or spline given by. the parametric curve obtained by applying the sine or Bezier or NURBS or spl ine interpolation method to the sampling sequence. respectively.
76 The device according to claim '72 which further includes the means for calculating the axis of the duct given by a vectorial function whose components are polynomials of the smallest degree with which an acceptable adjustment to the points of the axis sampling sequence is obtained.
77 The device according to any of claims 73 to 76 which further includes 101> differential and numerical calculation means for carrying out the analysis of the axis of the duct and comprising evaluating: a) the length of the duct, which is obtained by calculating the length of the curve that represents the cje adopted as a model and between the planes analogous to the conduit that delimit it in its two ends: oars; b) the orientation at each point of the axis and the average orientation at the points of the sampling sequence; e) instantaneous, average and maximum curvatures; d) instant torsions. average and maximum; e) the location of the points of maximum and minimum curvature and torsion; f) the Euclidean distances from said singular points to each end of the axis; and g) the distances along the axis between said singular points and each end of the axis.
The device according to any one of claims 73 to 76, which further includes the means for obtaining the image of the orthogonal cut to the axis of the conduit at an arbitrary point p thereof, together with the corresponding metadata, and the corresponding orthogonal binary mask. and that includes the following stages: a) calculate the tangent in p to the mathematical function adopted as the axis model; b} determine the equation of the plane perpendicular to the axis of the conduit enp; c) calculate the multiplanar reconstruction qu: e contains said plane perpendicular to the axis in p, together with the corresponding metadata, and obtain the orthogonal sectional image of the duct at point p along with its metadata. that it is equal to the image of the cut corresponding to said point in said multiplanar reconstruction, together with its metadata; d) apply to the multiplanar reconstruction the segmentation method of this invention and obtain the orthogonal binary mask of the duct in the cut perpendicular to its axis at point p that will be equal to the binary matrix corresponding to point p included in the obtained binary hyperpentry .
79 The device according to claims 65 or 78 which further includes the means for obtaining a plurality of 2D measurements of the section of the duct by an arbitrary point p of the axis thereof, each of which will be axial or orthogonal according to the section on the one obtained is axial or orthogonal respectively. measurements calculated all of them in the image coordinate system or in the physical coordinate system determined by the metadata. and which include one or more of the following measurements: a) sectional area of the duct in the cut by point p. obtained by adding the areas of the pixels marked as "conduit" in the corresponding binary mask of the conduit, the area of a pixel being equal to the product of its base by its height expressed in physical coordinates by the tomography metadata: b ) estimated area of the section of the canal in the orthogonal plane, obtained by multiplying the corresponding area in the axial section by the cosine of the angle between the axis of the tomography and that of the canal; e) maximum or minimum sectional diameter length. calculated as the maximum or minimum distance - respectively - between the outer edges of all the pairs that may form the pixels located on the periphery of the duct section; d) orientation of the duct section, given by the orientation of its maximum diameter: or e) moments of the section. standards. Central, centralized normalized and invariant. calculated according to universally accepted definitions; f) eccentricity of the section. calculated as the quotient between the distance between the foci and the length of the major axis of an ellipse with nominal second-order central moments equivalent to those of the duct section in the cut image. or by any other available procedure; g) elongation of the duct section, calculated by subtracting from the unit the quotient between the minor and major sides of the envelope rectangle of the section in the chosen coordinate system, or by any other available method; h) circularity of the s4 ~ ction. calculated as the inverse of 2p multiplied by the ratio between the square of the area and the sum of the variances in the two axes; i) section contour line, calculated as the set of points of contact with the bottom of the pixels belonging to the conduit: j) the plurality of indexes of the section that offers the geometric morphometry calculated directly on the contour line of the conduit or starting from marks or '"landmarks" derived from it.
80 The device according to claims 77 or 79 which also includes the means for obtaining statistics of each variable measured according to said claims by: a) the repeated carrying out of the measurements of said variable at all points p, of the Sampling sequence of the axis, or other set of axis points that is of interest: b) the calculation of the statistics of the measurements obtained. which include one or more of the following indices: i) the indices of central tendency and variability, ii) their extreme, minimum and maximum values; iii) determine the points of the axis where said extreme values are located; and iv) obtain the curves of each variable that express the measured value in ordinates and in abscissa the axis of the tomography or the rectified curve of the axis of the conduel.
81 The device according to any one of claims 45 to 80 which further includes the means for a) calculating the corresponding contour line from each binary axial or orthogonal sectional mask. ~ Axial or orthogonal equation - respectively ~ in the coordinate system chosen, line that as such has a null thickness and that is integrated by the set of contact points between the pixels of the duct and those of the bottom expressed in the physical coordinate system; b) generate from it the marks ("Iandmarks") and matrices of reference points in physical coordinates established in the metadata of the contour of the binary image of the conduit in the orthogonal cut to its axis, and represented in the corresponding image binary sectional; and e) calculate a plurality of geometric morphometry indexes 20
based on these brands.
82 The device according to any of claims 45 to 80 which further includes the
means to: a) obtain the marks ("Iandmarks") and matrices of reference points. in the physical coordinate system established in the metadata, based on solid 3-D models
or of the surface model or of the axis model of the duct; and b) calculate a plurality 3D geometric morphometry indices based on said marks and reference matrices.
The device according to any of claims 45 to 82. repeatedly applied to successive tomographs of the same conduit obtained over time, which also includes the means to a) calculate evolutionary indices of each variable and b) generate the curves and other formal models representative of the evolution of these variables over time.
The device according to claim 78 which further includes the means for creating the rectified tomography. following these stages; a) obtain the image of the oltogonal section to the axis of the duct at all points in the sampling sequence of the axis; e) position each image thus obtained in a plane perpendicular to the axial axis that cuts through the z coordinate
15 corresponding to that image and moving it in the indicated plane so that the centroid of the duct section coincides with the axis, thus obtaining the rectified tomography of the duct; d) record it in permanent memory: e) represent it graphically for visual analysis and metric morpho.
The device according to any one of claims 45 to 84. which further comprises the means for local or remote mode a) pennanently registering the models obtained; or b) make them accessible to the user, visualizing the models in virtual reality devices. screens or paper. in its entirety or by cuts: isolated or integrated into a volumetric rendering or other representation of its environment: in video mode. static figures, graphics. OR
25 interactive figures allowing the user to select the observation point, the colors and the nature and position of the light sources: oc) print in 20 or 3D of the results chosen by the user: od) generate digital copies on storage media Fixed or portable.
86 The device according to any of claims 45 to 85 further comprising the
30 means for local or remote mode a) record the graphic and numerical results obtained; or b) make them accessible to the user, displaying them on screens, paper or other physical or volatile physical media. in its entirety or in parts; In video mode, static figures. graphics. or interactive figures; or c) print the results chosen by the user; or d) generate digilalt copies: s in fixed or deliverable storage media.
87 The device according to any of claims 45 to 86 which further includes the means for modifying the fontate of the models obtained in such a way that they can be processed by 3D printers, numerical control machines or other means of machining or producing objects from numerical models capable of generating the surface model
5 a custom implant made in the type of material indicated by an expert or a physical representation of the duct model in a plurality of scales and materials.
88 The device according to any of claims 45 to 87, which further includes the means for calculating trajectories for guiding robotic devices, including surgical ones, 10 in cleaning, draining, filling, repair or surgery of the ducts represented in the surface model.
89 A system comprising the device according to any one of claims 45 to 88, and a wcb server. connected by a network to said device and also to the internet or other telematic networks; said server has the hardware and server and communications programs that allow it to) serve web pages, accept remote user connections intemet, receive the sending by the remote user of a request for analysis and modeling and receive the volume tomographic sent by the remote user containing the conduit OR conduits to be analyzed; b) transfer the request and the tomographic volume to that quoted
20 device for the analysis of the dimensions and morphometry of the duct; c) receive the results in electronic format from the device; d) transfer them to the applicant through the telematic channel or another that he has selected (: activated; and e) inform of the completion of the process other computer systems.

类似技术:

---

公开号 | 公开日 | 专利标题

---

ES2595210T3|2016-12-28|Method and apparatus for processing medical images

---

Echegaray et al.2015|Core samples for radiomics features that are insensitive to tumor segmentation: method and pilot study using CT images of hepatocellular carcinoma

---

Masood et al.2019|Automatic choroid layer segmentation from optical coherence tomography images using deep learning

---

Chen et al.2009|Automated ventricular systems segmentation in brain CT images by combining low-level segmentation and high-level template matching

---

CN106659424A|2017-05-10|Medical image display processing method, medical image display processing device, and program

---

Schreibmann et al.2014|Multiatlas segmentation of thoracic and abdominal anatomy with level set‐based local search

---

Ahammad et al.2018|Image processing based segmentation techniques for spinal cord in MRI

---

US9984462B2|2018-05-29|Disease characterization from fused pathology and radiology data

---

Wallner et al.2018|Clinical evaluation of semi-automatic open-source algorithmic software segmentation of the mandibular bone: Practical feasibility and assessment of a new course of action

---

Peng et al.2014|A region‐appearance‐based adaptive variational model for 3D liver segmentation

---

Poh et al.2012|Automatic segmentation of ventricular cerebrospinal fluid from ischemic stroke CT images

---

Cao et al.2016|An automatic breast cancer grading method in histopathological images based on pixel-, object-, and semantic-level features

---

Zheng et al.2016|Feature learning based random walk for liver segmentation

---

ES2608037A1|2017-04-05|System and method for segmentation and automated analysis of the three-dimensional structure of ducts in computed tomography images |

---

WO2021238438A1|2021-12-02|Tumor image processing method and apparatus, electronic device, and storage medium

---

Jimenez-Pastor et al.2020|Automated vertebrae localization and identification by decision forests and image-based refinement on real-world CT data

---

Cadotte et al.2015|Spinal cord segmentation by one dimensional normalized template matching: a novel, quantitative technique to analyze advanced magnetic resonance imaging data

---

Lenga et al.2018|Deep learning based rib centerline extraction and labeling

---

Memiş et al.2021|A novel approach for computerized quantitative image analysis of proximal femur bone shape deformities based on the hip joint symmetry

---

Zheng et al.2013|Adaptive segmentation of vertebral bodies from sagittal MR images based on local spatial information and Gaussian weighted chi-square distance

---

Serrurier et al.2012|Robust femur condyle disambiguation on biplanar X-rays

---

Babu et al.2019|A review on acute/sub-acute ischemic stroke lesion segmentation and registration challenges

---

Roy et al.2015|A useful approach towards 3D representation of brain abnormality from its 2D MRI slides with a volumetric exclamation

---

Jin et al.2019|How many models/atlases are needed as priors for capturing anatomic population variations?

---

Deng et al.2017|Semi-automatic segmentation of paranasal sinuses from CT images using active contour with group similarity constraints

同族专利:

---

公开号 | 公开日

---

ES2608037B1|2018-01-26|

引用文献:

---

公开号 | 申请日 | 公开日 | 申请人 | 专利标题

---

WO2003058553A2|2001-12-27|2003-07-17|The Government Of The United States Of America, As Represented By The Secretary Of The Department Of Health And Human Services|Automated centerline detection algorithm for colon-like 3d surfaces|

---

WO2011135103A1|2010-04-30|2011-11-03|Katholieke Universiteit Leuven, K.U. Leuven R&D|Intracoronary optical coherence tomography images analysis|

---

US20130121549A1|2010-07-30|2013-05-16|Koninklijke Philips Electronics N.V.|Organ-specific enhancement filter for robust segmentation of medical images|

---

US20140330115A1|2011-07-21|2014-11-06|Carrestream Health, Inc.|System for paranasal sinus and nasal cavity analysis|CN108765415A|2018-06-14|2018-11-06|南京鼎瑞医疗器械有限公司|There is one kind shade management to monitor system|

---

法律状态:
2018-01-26| FG2A| Definitive protection|Ref document number: 2608037 Country of ref document: ES Kind code of ref document: B1 Effective date: 20180126 |

---

2021-09-28| FD2A| Announcement of lapse in spain|Effective date: 20210928 |

优先权:

---

申请号 | 申请日 | 专利标题

---

ES201500707A|ES2608037B1|2015-10-01|2015-10-01|System and method for automated segmentation and analysis of the three-dimensional structure of ducts in computed tomography images|ES201500707A| ES2608037B1|2015-10-01|2015-10-01|System and method for automated segmentation and analysis of the three-dimensional structure of ducts in computed tomography images|