Method, system and computer program for determining position and/or orientation parameters of an anatomical structure

Abstract

Disclosed is a computer-implemented method of determining one or more position and/or orientation parameters of an anatomical structure of a body portion. The anatomical structure has a longitudinal shape defining a longitudinal axis. The method includes generating and/or reading, by a data processing system, volumetric data of at least a portion of a subject. The method further includes generating and/or reading, by the data processing system, a deformable template which provides an estimate for a location of the longitudinal axis in the portion of the subject. The method further includes matching, by the data processing system, the deformable template to the volumetric data, thereby obtaining a matched template. The matching comprises using one or more internal energy functions and one or more external energy functions for optimizing an objective function. The method further includes determining, by the data processing system, the at least one position and/or orientation parameter based on the matched template.

Claims

1. A computer-implemented method executed on one or more processors for determining one or more position and/or orientation parameters of an anatomical structure of a body portion, wherein the anatomical structure has a longitudinal shape defining a longitudinal axis, comprising: generating and/or reading, by the one or more processors, volumetric data of at least a portion of a subject; generating, by the one or more processors, a deformable template which provides an estimate for a location of the longitudinal axis in the portion of the subject; matching, by the one or more processors, the deformable template to the volumetric data, thereby obtaining a matched template; wherein the matching comprises using one or more internal energy functions and one or more external energy functions for optimizing an objective function; and determining, by the one or more processors, the one or more position and/or orientation parameters based on the matched template; wherein the deformable template represents a plurality of curves, each of which being an approximation for the longitudinal axis, and wherein each of the curves is at least partially defined using a plurality of position adjustable control points; wherein the deformable template comprises, for each of the control points, a discrete search space having a plurality of search space points for positioning the adjustable control points; wherein: (a) for at least a portion of the control points, the search spaces of the control points are mutually non-overlapping; and/or (b) for one or more of the control points, the search space points for the respective control point form a plane or curved surface.

2. The method of claim 1, further comprising: registering an atlas with the volumetric data; and determining a prototype of the deformable template based on the registered atlas.

3. The method of claim 1 further comprising determining, for each of a plurality of voxels, one or more values of an orientation probability density function of an orientation of the longitudinal axis at the respective voxel.

4. The method of claim 1 wherein the anatomical structure is at least a portion of a nerve fiber or at least a portion of a bundle of nerve fibers.

5. The method of claim 1 wherein at least one of the one or more internal energy functions depends on a distance of one or more control points of the deformable template from a prototype of the deformable template.

6. The method of claim 1 wherein at least one of the one or more external energy functions is a function of the locations of the control points and further a function of radius values of the anatomical structure at the locations.

7. The method of claim 1 wherein the objective function comprises a sum of subfunctions, which depends on all control points, wherein each of the subfunctions: depends on positions of directly consecutive control points, as seen along the curve; and is independent from the remaining control points.

8. The method of claim 7 further comprising: determining, for each search space point of a first one of the control points, an optimized energy value based on sections of the curves, which end at the respective search space point; wherein the optimized energy value is determined using one or a sum of the subfunctions, which corresponds to the curve sections.

9. A computer-implemented method executed on one or more processors for determining one or more position and/or orientation parameters of an anatomical structure of a body portion, wherein the anatomical structure has a longitudinal shape defining a longitudinal axis, comprising: generating and/or reading, by the one or more processors, volumetric data of at least a portion of a subject; generating, by the one or more processors, a deformable template which provides an estimate for a location of the longitudinal axis in the portion of the subject; matching, by the one or more processors, the deformable template to the volumetric data, thereby obtaining a matched template; wherein the matching comprises using one or more internal energy functions and one or more external energy functions for optimizing an objective function; and determining, by the one or more processors, the one or more position and/or orientation parameters based on the matched template; wherein: (a) at least one of the one or more internal energy functions depends on a distance of one or more control points of the deformable template from a prototype of the deformable template; and/or (b) the deformable template represents a plurality of curves, each of which being an approximation for the longitudinal axis, and wherein each of the curves is at least partially defined using a plurality of position adjustable control points; and at least one of the one or more external energy functions is a function of the locations of the control points and further a function of radius values of the anatomical structure at the locations.

10. The method of claim 9, further comprising: registering an atlas with the volumetric data; and determining a prototype of the deformable template based on the registered atlas.

11. The method of claim 9, further comprising determining, for each of a plurality of voxels, one or more values of an orientation probability density function of an orientation of the longitudinal axis at the respective voxel.

12. The method of claim 9, wherein the anatomical structure is at least a portion of a nerve fiber or at least a portion of a bundle of nerve fibers.

13. A computer-implemented method executed on one or more processors for determining one or more position and/or orientation parameters of an anatomical structure of a body portion, wherein the anatomical structure has a longitudinal shape defining a longitudinal axis, comprising: generating and/or reading, by the one or more processors, volumetric data of at least a portion of a subject; generating, by the one or more processors, a deformable template which provides an estimate for a location of the longitudinal axis in the portion of the subject; matching, by the one or more processors, the deformable template to the volumetric data, thereby obtaining a matched template; wherein the matching comprises using one or more internal energy functions and one or more external energy functions for optimizing an objective function; and determining, by the one or more processors, the one or more position and/or orientation parameter based on the matched template; wherein the deformable template represents a plurality of curves, each of which being an approximation for the longitudinal axis, and wherein each of the curves is at least partially defined using a plurality of position adjustable control points; wherein the objective function comprises a sum of subfunctions, which depends on all control points, wherein each of the subfunctions: depends on positions of directly consecutive control points, as seen along the curve; and is independent from the remaining control points; and further comprising: determining, for each search space point of a first one of the control points, an optimized energy value based on sections of the curves, which end at the respective search space point; wherein the optimized energy value is determined using one or a sum of the subfunctions, which corresponds to the curve sections.

14. The method of claim 13, further comprising: determining, for each search space point of a second one of the control points, an optimized energy value based on sections of the curves, which end at the respective search space point and which pass through the search space of the first control point; wherein the second control point is a directly consecutive control point to the first control point; and wherein the optimized energy value is determined based on the optimized energy values determined for the first control point and further based on the subfunction which depends on the first and second control point.

15. The method of claim 13, further comprising: registering an atlas with the volumetric data; and determining a prototype of the deformable template based on the registered atlas.

16. The method of claim 13, wherein the anatomical structure is at least a portion of a nerve fiber or at least a portion of a bundle of nerve fibers.

17. A non-transitory computer readable media comprising instructions for determining one or more position and/or orientation parameters of an anatomical structure of a body portion, wherein the anatomical structure has a longitudinal shape defining a longitudinal axis, which when executed by at least one processor, causes the at least one processor to: generate and/or reading, by the at least one processor, volumetric data of at least a portion of a subject; generate, by the at least one processor, a deformable template which provides an estimate for a location of the longitudinal axis in the portion of the subject; match, by the at least one processor, the deformable template to the volumetric data, thereby obtaining a matched template; wherein the matching comprises using one or more internal energy functions and one or more external energy functions for optimizing an objective function; and determine, by the at least one processor, the one or more position and/or orientation parameters based on the matched template; wherein the deformable template represents a plurality of curves, each of which being an approximation for the longitudinal axis, and wherein each of the curves is at least partially defined using a plurality of position adjustable control points; wherein the deformable template comprises, for each of the control points, a discrete search space having a plurality of search space points for positioning the adjustable control points; wherein: (a) for at least a portion of the control points, the search spaces of the control points are mutually non-overlapping; and/or (b) for one or more of the control points, the search space points for the respective control point form a plane or curved surface.

18. A system for determining one or more position and/or orientation parameters of an anatomical structure of a body portion, wherein the anatomical structure has a longitudinal shape defining a longitudinal axis, comprising: a radiation treatment apparatus comprising a treatment beam source and a patient support unit; at least one processor executing instructions to: generate and/or read, by the at least one processor, volumetric data of at least a portion of a subject; generate, by the at least one processor, a deformable template which provides an estimate for a location of the longitudinal axis in the portion of the subject; match, by the at least one processor, the deformable template to the volumetric data, thereby obtaining a matched template; wherein the matching comprises using one or more internal energy functions and one or more external energy functions for optimizing an objective function; and determine, by the at least one processor, the one or more position and/or orientation parameter based on the matched template; wherein the deformable template represents a plurality of curves, each of which being an approximation for the longitudinal axis, and wherein each of the curves is at least partially defined using a plurality of position adjustable control points; wherein the deformable template comprises, for each of the control points, a discrete search space having a plurality of search space points for positioning the adjustable control points; wherein: (a) for at least a portion of the control points, the search spaces of the control points are mutually non-overlapping; and/or (b) for one or more of the control points, the search space points for the respective control point form a plane or curved surface; wherein the at least one processor is operably coupled to the radiation treatment apparatus for issuing a control signal to the radiation treatment apparatus for controlling, on the basis of the determined position and/or orientation parameter of the anatomical structure, at least one of the operation of the treatment beam source or the position of the patient support unit.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) In the following, the present disclosure is described with reference to the appended Figures which give background explanations and represent exemplary embodiments of the invention. The scope of the invention is however not limited to the specific features disclosed in the context of the Figures, wherein

(2) FIG. 1 is a schematic illustration of a medical system according to an exemplary embodiment;

(3) FIG. 2 is a flow chart illustrating an exemplary method for determining one or more position and/or orientation parameters of an anatomical structure of a body portion performed by the medical system according to the exemplary embodiment shown in FIG. 1;

(4) FIG. 3A is a cross-section through volumetric data which are acquired and analyzed using the medical system according to the exemplary embodiment, which is shown in FIG. 1;

(5) FIG. 3B is a schematic illustration showing an atlas which is registered to the volumetric data according to the exemplary method, which is illustrated in FIG. 2;

(6) FIG. 4 is a schematic illustration of the determination of an estimate for the longitudinal axis in the exemplary method, which is illustrated in FIG. 2;

(7) FIGS. 5A and 5B schematically illustrate a deformable template which is used in the exemplary method, which is illustrated in FIG. 2;

(8) FIGS. 6A and 6B schematically illustrate different curves generated using the deformable template in the exemplary method, which is illustrated in FIG. 2;

(9) FIG. 7 schematically illustrates a matched template and an estimate for an extent of the anatomical structure, which are determined in the exemplary method, which is illustrated in FIG. 2;

(10) FIG. 8 schematically illustrates a comparison between the volumetric data and the determined estimate for the extent of the anatomical structure which is determined using the exemplary method, which is illustrated in FIG. 2; and

(11) FIG. 9 schematically illustrates a second example of a deformable template which is used by the imaging system in a second exemplary method; and

(12) FIG. 10 schematically illustrates a third example of a deformable template which is used by the data processing system in a third exemplary method.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

(13) In the following, a short description of the specific features of the present disclosure is given which shall not be understood to limit the invention only to the features or a combination of the features described in this section.

(14) FIG. 1 is a schematic illustration of a medical system 1 according to an exemplary embodiment. The medical system 1 includes a data processing system 2 and a data acquisition system 7, which is in signal communication with the data processing system 2. In the exemplary embodiment shown in FIG. 1, the data acquisition system 7 is configured as a magnetic resonance tomography system. It is to be understood, however, that the present disclosure is not limited to data acquisition systems, which include magnetic resonance tomography systems, but encompasses all embodiments of data acquisition systems, which are configured to acquired volumetric data from a subject under inspection.

(15) The data processing system 2 may be a stand-alone computer or may be configured as a distributed computer system which uses a computer network 3, such as the Internet or a local area network (LAN). The data processing system 2 includes a display device 4 and input devices, such as a keyboard 5 and a computer mouse 6 allowing user interaction via a graphical user interface of the data processing system 2.

(16) FIG. 2 is a flow chart illustrating an exemplary method 100 performed using the medical system 1 which is schematically illustrated in FIG. 1.

(17) The data processing system 2 is configured to read and/or generate S101 volumetric data of at least a portion of the subject. By way of example, the data processing system 2 may read the volumetric data from the data acquisition system 7 and/or may be configured to generate the volumetric data from data from the data acquisition system 7 and/or from signals received from the data acquisition system 7.

(18) In the exemplary embodiment, which is shown in FIG. 1, the data acquisition system 7 is configured to perform diffusion tensor magnetic resonance imaging, which is a technique of diffusion-weighted magnetic resonance imaging and which has been shown to be able to depict brain white matter bundles.

(19) In order to perform diffusion tensor magnetic resonance imaging, the data processing system 2 or the data acquisition system 7 (such as a further data processing system which is part of the data acquisition system 7) is configured to generate, depending on data acquired from a subject under inspection, for each of a plurality of voxels, values of a diffusion tensor. The diffusion tensor D is a [3×3] symmetric matrix:

(20) $D=\left[\begin{array}{ccc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}& D_{\mathrm{xz}}\\ D_{\mathrm{yx}}& D_{\mathrm{yy}}& D_{\mathrm{yz}}\\ D_{\mathrm{zx}}& D_{\mathrm{zy}}& D_{\mathrm{zz}}\end{array}\right],$

(21) where each of the three diagonal elements (D.sub.xx, D.sub.yy, D.sub.zz) represents an apparent diffusion coefficient measured along one of the laboratory axes (i.e. the x-, y- and the z-axis). Each of the six off-diagonal terms (D.sub.xy, D.sub.yz, etc.) reflect the correlation of random motion between a pair of directions, which correspond to two different laboratory axes.

(22) For the special case of perfect isotropic diffusion (such as seen in pure liquids), the off-diagonal elements are all zero and the diagonal elements are all the same and equal the single diffusion coefficient, D.sub.0, for the isotropic material (i.e. D.sub.xx=D.sub.yy=D.sub.zz=D.sub.0). For anisotropic diffusion, however, as it occurs in diffusion tensor magnetic resonance imaging of brain white matter, for at least a portion of the voxels, the diagonal elements of the diffusion tensor are unequal and the off-diagonal elements cannot be ignored.

(23) FIG. 3A schematically illustrates a two-dimensional cross-section 10 through a volumetric data set generated using the medical system 1 (shown in FIG. 1). In FIG. 3A, for each of the voxels, the grayscale value of the two-dimensional cross-section 10 indicates a level of fractional anisotropy (commonly abbreviated as FA) at the respective voxel, which is determined using values of the diffusion tensor at the respective voxel.

(24) The FA level is a scalar value between 0 and 1 that describes the degree of anisotropy of a diffusion process. A value of 1 (which represents a theoretical extreme) means that diffusion occurs only along one axis and is fully restricted along all other directions. The FA level is calculated from the eigenvalues λ.sub.1, λ.sub.2 and λ.sub.3 of the diffusion tensor using the equation:

(25) $\mathrm{FA}=\sqrt{\frac{3}{2}}\frac{\sqrt{{\left({\lambda }_1-\hat{\lambda }\right)}^2+{\left({\lambda }_2-\hat{\lambda }\right)}^2+{\left({\lambda }_3-\hat{\lambda }\right)}^2}}{\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}},$

(26) with {circumflex over (λ)}=(λ.sub.1+λ.sub.2+λ.sub.3)/3 being the mean value of the eigenvalues.

(27) Therefore, a high FA level means a more unidirectional flow whereas a low FA level means equal water movement in all directions. As such, as a general rule, image regions of brighter grayscale values indicate greater white matter anisotropy which may be caused by highly organized fiber bundles which are present at the image region.

(28) It is to be noted that in the exemplary method, the FA values are used for illustrating the performance of the invention so that is not necessary to determine values of the FA for performing the exemplary method. On the other hand, it is also conceivable that one or more of the external energy functions of the deformable template are determined using values of FA.

(29) It is known in the prior art to use deterministic or probabilistic fiber tracking algorithms to determine the extent of brain white matter fiber bundles based on diffusion tensor imaging data. These methods, however, have shown to suffer from significant limitations.

(30) Specifically, deterministic fiber tacking methods do not account for uncertainty in the determined pathway so that valuable information about the underlying diffusion tensor imaging data is lost. Further, deterministic fiber tracking methods do not consider pathways, which pass through small regions that violate deterministic rules and which may be caused by measurement artefacts, noise and unresolved features, such as crossing fiber tracts.

(31) While probabilistic tractography algorithms expand the pathway search space beyond that explored by deterministic algorithms and explicitly represent uncertainty in the data, probabilistic tractography algorithms also suffer from the limitation that they do not yield an accurate probability of brain connections. Therefore, like deterministic tracking algorithms, prior knowledge of anatomy fiber tracts is important for distinguishing between fiber tracts of interest and tracks that follow improbable routes or suggest non-existent connections between brain areas. Moreover, probabilistic tractography algorithms typically require intensive computation which inhibits its application in routine clinical tasks.

(32) However, as is explained in the following, in view of the above-described limitations of conventional fiber tracking algorithms, the present inventors have recognized that it is possible to provide an improved method for fiber tracking. As is also explained further below, the inventors also have recognized that it is possible to provide an improved analysis procedure for other anatomical structures than brain white matter fiber bundles, such as vascular anatomical structures, and particularly blood vessels.

(33) As is indicated in the flow chart of FIG. 2, the data processing system is configured to read and/or to generate S102 data representing an estimate for a longitudinal axis of the anatomical structure.

(34) By way of example, the data processing system is configured to use segmentation data, which was read and/or generated by the data processing system and to determine the estimate for the longitudinal axis depending on the segmentation data. The segmentation data indicates, for each of a plurality of voxels of the volumetric data, an estimate whether the respective voxel represents a part of the anatomical structure. Therefore, the segmentation data may include a voxel object. In FIG. 3B, the segmentation data are schematically indicated by a white contour 11 within the two-dimensional cross section 10 through the volumetric data.

(35) The segmentation data may be generated by the data processing system or by an external computer system. The segmentation data may be generated using an atlas-based segmentation procedure. However, it is also conceivable that, additionally or alternatively, other segmentation procedures than atlas-based segmentation procedures are used. By way of example, it is conceivable that the segmentation data are generated using an artificial neural network. The segmentation procedure may be automatic or semi-automatic (i.e. requiring user intervention).

(36) The atlas may be generated by integrating segmentation data from multiple segmented images. The multiple images may be generated from different individuals. The atlas may thereby represent an average shape atlas.

(37) The atlas may be registered with at least a portion of diffusion-weighted magnetic resonance data and/or with at least a portion of non-diffusion-weighted magnetic resonance data.

(38) By way of example, the atlas may be registered (in particular directly registered) with magnetic resonance data using a b0 scan (i.e. a non-diffusion weighted scan). The b0 scan may be part of a diffusion tensor imaging data set. The volumetric image data may include the diffusion tensor imaging data set or the volumetric image data may be generated based on the diffusion tensor imaging data set.

(39) Alternatively, it is also conceivable that the atlas is indirectly registered with the diffusion-weighted magnetic resonance data. By way of example, the atlas may be registered to anatomical magnetic resonance data which in turn is registered with the diffusion-weighted magnetic resonance data.

(40) Examples for anatomical magnetic resonance data are but are not limited to: T1 and T2-weighted magnetic resonance data and Constructive interference in Steady State (CISS) data.

(41) In FIG. 3B, the segmentation data are schematically illustrated as a white contour line 12 in the cross-section 10 through the volumetric data. As can be seen from FIG. 3B, the extent of the anatomical structure, which is derived using atlas registration, is significantly different from the extent of the anatomical structure, as suggested by the levels of fractional anisotropy. However, it has been shown by the inventors that the procedure, which is described in detail in the following, generates an estimate for the extent of the anatomical structure, which is more consistent with the diffusion tensor imaging data.

(42) As is schematically illustrated in FIG. 4, the data processing system may be configured to generate, depending on the segmentation data 11, the estimate 12 for the longitudinal axis.

(43) By way of example, determination of the estimate for the longitudinal axis may include successively removing outer layers of the voxel object until a line-shaped object is obtained having a diameter of less than five voxels, or substantially one voxel. Additionally or alternatively, the determination of the estimate for the longitudinal axis may include determining a spline curve representing the estimate for the longitudinal axis. However, it is noted that the present disclosure is not limited to one or a combination of these procedures of determining the longitudinal axis.

(44) As is indicated in FIG. 2, the data processing system is configured to generate S103 a deformable template. Specifically, as is illustrated in FIG. 5A, the data processing system is configured either to use the estimate for the longitudinal axis 12 as a prototype 14 for the deformable template 13 or to generate the prototype 14 based on the estimate for the longitudinal axis 12 (e.g. by approximating the estimate as a polygonal chain). The deformable template 13 is configured as a parametric deformable template representing a set of deformed curves, which are uniquely described by values of a set of parameters. In other words, the geometrical shape, position and/or orientation of the curves of the deformable template 13 can be changed by using different parameter values. The prototype 14 describes only one of the plurality of curves represented by the deformable template 13.

(45) FIG. 5B schematically illustrates the deformable template 13 in more detail. The deformable template 13 includes a plurality of position adjustable control points 15a, . . . 15i, (shown as circles in FIG. 5B), which are connected by line segments so that each of the line segments connects two of the position adjustable control points 15a, . . . 15i. The line segments which are shown in FIG. 5B together form the prototype 14. By way of example, the data processing system determines the locations of the position adjustable control points 15a, . . . 15i of the prototype 14 so that a distance between neighboring position adjustable control points is substantially constant.

(46) Each of the position adjustable control points, 15a, . . . 15i, is adjustable within a discrete search space 16a, . . . 16i, having a plurality of search space points (indicated in FIG. 5B as black dots) so that separate search spaces are provide for each of the position adjustable control points 15a, . . . 15i. In the example, which is shown in FIG. 5B, the search spaces are mutually non-overlapping. Hover, it is also conceivable that at least a portion of the search spaces are overlapping, such as in the exemplary embodiment, which is described further below in connection with FIG. 9.

(47) Therefore, the positions of the position adjustable control points within their respective search spaces represent the set of parameters, which defines the shape of the curves represented by the deformable template. Thereby, different positions of the control points within their respective search spaces generate different curves, as is schematically illustrated in FIGS. 6A and 6B.

(48) However, it is noted that the present disclosure is not limited to deformable templates having a polygonal chain. By way of example, the curve may be a spline curve having a degree greater than 1. A spline curve of degree 1 is a polygonal chain. Additionally or alternatively, it is also conceivable that the deformable template having a Bezier curve and/or a NURBS (non-uniform rational B-spline) curve.

(49) For the sake of easy understanding, in the schematic illustration of FIG. 5B, the search space points, which are part of a common search space, are connected by a line. In the illustrated exemplary embodiment, the deformable template is configured so that each of the discrete search spaces 16a, . . . 16i forms a plane surface, which is oriented perpendicular to a tangent to the prototype 14 at a location of one of the position adjustable control points of the prototype 14. By way of example, the surfaces which are spanned by the search space points of a discrete search space may be in the form of a disk. It is, however, conceivable that the discrete search spaces have different shapes and/or have a non-planar geometry.

(50) As is further indicated in the flowchart of FIG. 2, the processing system is configured to match S104 the deformable template to the volumetric data in order to obtain a matched template. The matched template is a curve of the deformable template, which optimizes (i.e. maximizes or minimizes) an objective function. The inventors have shown that through such a matching process, a curve can be identified, which represents an estimate for the location of the longitudinal axis of the anatomical structure.

(51) FIGS. 6A and 6B show two curves of the deformable template, which are generated using different sets of parameter values. The curve 18, which is shown in FIG. 6A has a higher degree of smoothness compared to the curve 23, which is shown in FIG. 6B. Further, as measured within the search space of each of the control points, a predominant portion of the control points of the curve 18 of FIG. 6 have a smaller distance from the prototype 14 than the corresponding control points of the deformed template of FIG. 6B. It has been shown that a more reliable estimate for the longitudinal axis of the anatomical structure can be obtained if the matching process favors deformed templates having a smooth appearance and which are located at a small distance from the prototype 14. As will be discussed in the following, this is achieved by the objective function of the deformable template, which is optimized during the matching process.

(52) In the exemplary embodiment, the matching process includes maximizing an objective function, which is the sum of an internal energy function and an external energy function:
OF=E.sub.int+E.sub.ext

(53) For each of the curves generated by the deformable template, the internal energy function depends on one or more intrinsic geometric parameters of the respective curve and further on one or more parameters indicative of a distance of the curve from the prototype. The internal energy function is independent from the volumetric data.

(54) In the exemplary embodiment, the internal energy function E.sub.int is defined as:
E.sub.int=E.sub.dist+E.sub.bend,

(55) with E.sub.dist being the distance energy function, which penalizes large distances from the prototype and E.sub.bend being the bending energy, which penalizes sharp bends.

(56) By way of example, the distance energy E.sub.dist, may be expressed as:

(57) $E_{\mathrm{dist}}=a_{\mathrm{int}}\log \underset{i=1}{\overset{N}{.Math.}}\exp \left(\frac{d_i^2}{{\sigma }^2}\right)=\underset{i=1}{\overset{N}{.Math.}}\frac{d_i^2}{{\sigma }^2};$

(58) with N being the number of control points, d.sub.i being the distance of the i-th control point from the corresponding control point position of the prototype. σ is a parameter, the value of which determines the degree to which small distances are penalized. The factor α.sub.int is a weighting factor that weights the relative contribution of the distance energy E.sub.dist to the internal energy E.sub.int.

(59) Further by way of example, the bending energy E.sub.bend may be expressed as:
E.sub.bend=α.sub.bendΣ.sub.i=1.sup.N−1f(θ.sub.i)

(60) wherein θ.sub.i is an angle of the i-th segment (which connects the i-th control point with the i+1-th control point) relative to the plane-shaped discrete search space of the i-control point. In an alternative embodiment, the angle θ.sub.i is the angle between the i-th segment and the i+1-th segment. f(θ.sub.i) is a function, which penalizes large angles. By way of example, if θ.sub.i is larger than 35 degrees, the function is not admitted as a candidate for the matched template (e.g. by setting f(θ.sub.i) to minus infinity). If θ.sub.i is smaller than 35 degrees, the curve is not penalized (e.g. by setting f(θ.sub.i) to a constant value of 0). The factor α.sub.bend is a weighting factor that weights the relative contribution of the bending energy E.sub.bend to the internal energy E.sub.int.

(61) As such, the internal energy E.sub.int penalizes a low degree of smoothness as well as large distances from the prototype.

(62) It is to be noted that the present disclosure is not limited to the above expressions for E.sub.int, E.sub.dist and E.sub.bend. Various modifications are conceivable regarding the functional form of the internal energy E.sub.int.

(63) The external energy function, on the other hand, depends on values of the volumetric data. The external energy function depends on an orientation probability density function of an orientation of the longitudinal axis of the anatomical structure. The orientation probability density function is determined for each of a plurality of voxels based on the diffusion tensor data at the respective voxel. At each of the position adjustable control points of a curve, the orientation probability density function is evaluated using one of the adjacent line segments. Specifically, the external energy E.sub.ext is calculated as:
E.sub.ext=logπ.sub.i=1.sup.N−1p(D.sub.i|t.sub.i)=Σ.sub.i=1.sup.N−1 log(d.sub.i|t.sub.i),

(64) wherein t.sub.i represents a line segment of the curve which connects the i-th control point with the i+1-th control point, N is the number of control points and D.sub.i is the diffusion tensor at the ith control point. Therefore, t.sub.i represents a tangent to the curve at the i-th control point. The function p(D.sub.i|t.sub.i) is the orientation probability density function, which indicates, for the direction of t.sub.i at the position of the i-th control point, a probability value that the longitudinal axis of the anatomical structure is oriented along the direction. In exemplary embodiments, in which the curve is not a polygonal chain, t.sub.i represents a tangent to the respective curve at the i-th control point. Therefore, determining the external energy function for a curve of the deformable template may include evaluating, for each of a number of points, the orientation probability density function at the respective point. The orientation probability density function may be evaluated for a direction, which corresponds to a tangent of the curve at the respective point.

(65) It is to be noted that the present disclosure is not limited to a specific functional form of the orientation probability density function. An example for an orientation probability density function is disclosed in the article “Modelling noise-induced fiber orientation error in diffusion-tensor MRI”, written by Philip Cook et al. and published in “2004 2nd IEEE International Symposium on Biomedical Imaging: Nano to Macro (IEEE Cat No. 04EX821)”, which is incorporated by reference herein in its entirety. This article introduces an orientation probability density function using the Watson Distribution on a sphere. A further example for an orientation probability density function is disclosed in the article “ConTrack: Finding the most likely pathways between brain regions using diffusion tractography”, written by Anthony J. Sherbondy et al. and published in the “Journal of Vision”(2008) 8 (9): 15, pages 1 to 16, which is incorporated by reference herein in its entirety. This article proposes using a Bingham distribution for determining the orientation probability density function.

(66) In exemplary embodiments, in which the volumetric data are generated using diffusion spectrum imaging (DSI), the determination of the orientation probability density function may include determining an inverse 3D-Fourier transformation based on the diffusion-weighted volumetric images acquired at the sampled q-space points. Further, in exemplary embodiments, in which the volumetric data are generated using single shell or multi-shell high angular resolution imaging (HARDI), the determining of the orientation probability density function may include using a Funk-Radon transform. One of such techniques is termed q-ball imaging and described in the article “Q-Ball Imaging”, written by David S. Tuch and published in “Magnetic Resonance in Medicine” 52:1358-1372 (2004), the contents of which is incorporated by reference herein in its entirety.

(67) In the exemplary embodiment, the curve, which maximizes the objective function is determined using an exhaustive search. A further, computationally efficient method for determining the matched template is described further below in connection with FIG. 9. FIG. 7 shows a matched template 19, which is a curve of the deformable template, which is obtained by optimizing the objective function. As can be seen from FIG. 7, the total energy function results in a matched template, which has no sharp bends. The matched template 19 represents an approximation for the longitudinal axis of the anatomical structure. Therefore, through the matching process, the data processing system has determined S105 (shown in FIG. 2) position and orientation parameters of the longitudinal axis of the anatomical structure.

(68) The data processing system uses the matched template 19 to determine an extent 20 of the anatomical structure, as it is schematically illustrated in FIG. 7. In order to determine the extent 20 of the anatomical structure, the data processing system uses prior knowledge about the diameter of the anatomical structure (e.g. the diameter of the cranial nerve). Specifically, the data processing system determines a surface, which surrounds the matched template 19 and which is located at a radius from the matched template 19, which corresponds to half of the known diameter of the anatomical structure. Furthermore, also a known length of the anatomical structure can be used to adapt the matched template 19. By way of example, one or both ends of the anatomical structure may be known. Specifically, it is known that all cranial nerves originate from the brain stem. If the shape of the search space of the first control point does not correspond to the shape of the brain stem, this can be corrected by adjusting the length of the matched template.

(69) FIG. 8 is a comparison between the determined extent 20 and the fractional anisotropy data in the cross-section 10 through the volumetric data. It can be seen from this comparison that the extent 20 which is determined based on the deformable template matches the extent as indicated by the fractional anisotropy levels to a high degree.

(70) FIG. 9 is a schematic illustration of a deformable template 21, which includes a comparatively high number of position adjustable control points. Due to the high number of discrete search spaces, the total number of curves defined by these search spaces is quite high. However, as is explained in the following, the inventors have shown that it is possible to determine the matched template 19 with a comparatively low computational effort.

(71) Specifically, the inventors have shown, that a computationally efficient determination of the matched template is possible using an objective function (OF), which is expressed as a sum of subfunctions ψ:

(72) $\mathrm{OF}\left(,\mathrm{.Math.},\right)=\underset{k=1}{\overset{N-1}{.Math.}}\psi \left(,\right),$

(73) with ν.sub.1, . . . ν.sub.N being the control point positions and N being the number of control points. Each of the subfunctions ψ depends on control point positions of a neighboring pair of the control points (i.e. a directly consecutive pair, as seen along the longitudinal axis), which are connected by a line segment of the curve and is independent from the positions of the remaining control points of the deformable template. It is noted that the objective function, which was discussed above in connection with the deformable template shown in FIG. 5 fulfills this requirement.

(74) As is explained in the following, the above form of the objective function

(75) OF allows determination of the matched template by successively optimizing partial objective functions POF.sub.j, for values of j=1 to N−1 (N being the number of control points):
POF.sub.j=Σ.sub.k=1.sup.jψ(,),

(76) wherein the optimization is processed in order from j=1 to N−1 and (shown in FIG. 9) is a given position of the first control point within its search space and may correspond to the first control point of the prototype 14.

(77) Specifically, for a curve, which starts at the given location of of the first control point, for each search space point in the search space of the second control point , the energy for the curve section between the first and the second control point can be expressed as:
POF.sub.1=ψ(,),

(78) i.e. the energy is determined using the subfunction, which corresponds to the curve section.

(79) As a next step, using the energy values assigned to the search space points within the search space of the second control point, the data processing system identifies, for each search space point in the discrete search space of the third control point, a search space point of the second control point, which maximizes the partial objective function POF.sub.2 using the given control point position of the first control point. Thereby, an optimized energy value and a corresponding optimized curve section can be assigned to each search space point in the discrete search space of the third control point.

(80) In other words, the optimized energy value is determined using a sum of the subfunctions (i.e. POF.sub.2), which correspond to the curve sections.

(81) Then, as a further step, using the energy values assigned to the search space points within the search space of the third control point, the data processing system identifies, for each search space point in the discrete search space of the fourth control point, a search space point within the search space of the third control point, which maximizes the partial objective function POF.sub.3 using the optimized energy values assigned to the search space points of the third control point. Thereby, an optimized energy value and a corresponding optimized section of a curve can be assigned to each search space point in the discrete search space of the fourth control point.

(82) In other words, for each of the search space points of the fourth control point, the optimized energy value is determined based on the optimized energy values determined for the third control point and further based on the subfunction ψ(,), which depends on the locations of the third and fourth control point.

(83) Using this procedure, it is possible to determine, for each search space point in the discrete search space of the last (i.e. N-th) control point, an optimized energy value and a curve, which starts at and which ends at the respective point in the discrete search space of the last control point.

(84) The data processing system then determines the point in the discrete search space of the last control point, which maximizes the objective function. The corresponding curve 19 (shown in FIG. 9) is then determined to be the matched template.

(85) As can further be seen from the deformable template which is illustrated in FIG. 9, a portion of the discrete search spaces overlap. The overlapping search spaces can lead to curves having a high bending energy, since the curve may change its direction by almost 180 degrees. The matching algorithm can be configured to rearrange for such curves the order of the search spaces in order to avoid sharp bends.

(86) FIG. 10 is a schematic illustration of a portion of a third example of a deformable template 13 which is used by the data processing system in third exemplary method for determining the longitudinal axis of a tubular anatomical structure, such as a vascular body (e.g. a blood vessel).

(87) As is explained in the following in detail, the deformable template 13 according to the third example is configured not only to determine an estimate for the longitudinal axis of the anatomical structure but also to determine, for each control point position of the matched template, an estimate for the radius of the anatomical structure.

(88) The volumetric data based on which the matched template is determined may be anatomical MRT data, such as T1 and/or T2 images. Additionally or alternatively, it is conceivable that one or a combination of the following data acquisition techniques are used: diffusion-weighted magnetic resonance imaging (e.g. using values of fractional anisotropy determined from values of the diffusion tensor), computed tomography, sonography, ultrasound imaging, and positron emission tomography. The volumetric data may represent a three-dimensional scalar field.

(89) In the method according to the third example, the internal energy is the same as mentioned above in connection with the method according to the first example.

(90) As is illustrated in FIG. 10, in the method according to the third example, the deformable template 13 is configured so that for each control point (such as the control points at locations and shown in FIG. 10), there are a plurality of search spaces provided (such as the search spaces 24a, 25a and 26a for the control point located at and the search spaces 24b, 25b and 26b for the control point located at ).

(91) For a given control point, each of the search spaces corresponds to a different radius value. By way of example, each of the search spaces 24a and 24b corresponds to a first radius value R.sub.1, each of the search spaces 25a and 25b corresponds to second radius value R.sub.2 and each of the search spaces 26a and 26b correspond to a third radius value R.sub.3.

(92) For the sake of easy illustration, in FIG. 10, the search spaces which are associated to a same control point are illustrated as being partially overlapping. However, in the deformable template 13 of the third example, at each location of a search point (such as the search point location which is shown in FIG. 10), there are points of three search spaces located, wherein each of the search spaces corresponds to a different radius value (R.sub.1, R.sub.2, R.sub.3).

(93) The deformable template 13 is configured to determine an estimate for an extent of the anatomical structure by providing an estimate for the longitudinal axis and estimates for radius values at a plurality of locations along the longitudinal axis. The radius values are limited to discrete values, such as the three values R.sub.1, R.sub.2, R.sub.3. It is noted in this regard that the number of radius values (and therefore the number of search spaces for each control point) may be greater than three or less than three. The number may be chosen depending on the anatomical structure, which is analyzed and/or depending on the required accuracy. The number of discrete radius values (and therefore the number of search spaces) as well as the values for the radii may be identical for each control point of the deformable template or may be different among the control points.

(94) In the method according to the third example, the data processing system determines for each of the search space points of a control point (i.e. for each of the search points in the plurality of search spaces which are associated with different radius values), a probability measure P which is a measure for the probability that the longitudinal axis passes through the respective search space point and has a radius value which is assigned to the respective search space point.

(95) The probability measure P may be determined using a Gradient Vector

(96) Flow vector field which is determined using the volumetric data and by solving a partial differential equation. The vectors of the Gradient Vector Flow vector field point to object edges of the anatomical structure. However, it is also conceivable that the values of the probability measure P are determined using the gradient vector field which is determined based on the volumetric data.

(97) Specifically, using the Gradient Vector Flow field, the probability measure P may be determined for the search space point located at a location as follows: an estimate for a cross-sectional plane is determined which passes through and which is approximately located perpendicular to the longitudinal axis. The estimate for the cross-sectional plane may be determined using a least squares fitting method based on the Gradient Vector Flow field. Additionally or alternatively, the estimate for the cross-sectional plane may be determined using the prototype for the deformable template. By way of example, the cross-sectional plane may be determined so that it passes through and is oriented perpendicular to a tangent to the prototype.

(98) Using the estimated cross-sectional plane, for each of the three search space points at the location (i.e. the search space point for the search space 24a, the search space point for the search space 24b and the search space point for the search space 24c), the probability measure P is determined using the flux of the Gradient Vector Flow field (or the gradient vector field) which is projected onto the cross-sectional plane and which passes through a circle having the radius of the respective search space point (i.e. the radius value R.sub.1, R.sub.2 or R.sub.3).

(99) The external energy E.sub.ext for a curve of the deformable template is then expressed using the determined values for the probability measure P. Specifically, a curve which is defined by the control points locations . . . and which has the radius values . . . at the respective control point locations ({circumflex over (R)}.sub.l being one of R.sub.1, R.sub.2 or R.sub.3), the external energy is expressed as:
E.sub.ext=Σ.sub.i=1.sup.NP(,{circumflex over (R)}.sub.l),

(100) with N being the number of control points. Therefore, the external energy function is a function of the locations of the control points and further a function of a radius value for each of the control points.

(101) Therefore, in a similar manner as has been described in connection with the method according to the second example, the objective function can be expressed as a sum of subfunctions:
OF(, . . . ,,, . . . )=Σ.sub.k=1.sup.N−1ψ(,),

(102) so that the process of matching the deformable template to the volumetric data can be carried out in a similar manner as has been described above in connection with the method according to the second example.

(103) In other words, each of the subfunctions depends on two radius values, each of which representing a radius of the anatomical structure at one of the consecutive control points.

(104) Specifically, as is explained in the following, the above form of the objective function OF allows determination of the matched template by successively optimizing partial objective functions POF.sub.j, for values of j=1 to N−1 (N being the number of control points):
POF.sub.j=Σ.sub.k=1.sup.jψ(,,,),

(105) wherein the optimization is processed in order from j=1 to N−1 and is a given position of the first control point within its search spaces and may correspond to the first control point of the prototype 14.

(106) For each search space point in the search spaces of the second control point , (i.e. for each location and radius value), the energy for the curve section between the first and the second control point can be expressed as:
POF.sub.1=ψ(,,,),

(107) wherein is the given position of the first control point within its search space and the Radius value can be obtained by optimizing the external energy function which corresponds to the curve section (i.e. E.sub.ext=P(,)+P(,)),

(108) As a next step, using the energy values assigned to the search space points of the second control point, the data processing system identifies, for each search space point of the third control point, a search space point (i.e. a location and a radius value) of the second control point, which maximizes the partial objective function POF.sub.2 using the given control point position and the given radius value of the first control point. Thereby, an optimized energy value and a corresponding optimized curve section (including the radius values at the control points of the optimized curve section) can be assigned to each search space point of the third control point.

(109) Then, as further step, using the energy values assigned to the search space points of the third control point, the data processing system identifies, for each search space point of the fourth control point, a search space point (i.e. a location and a radius value) of the third control point, which maximizes the partial objective function POF.sub.3 using the optimized energy values assigned to the search space points of the third control point. Thereby, an optimized energy value and a corresponding optimized curve section (including the radius values at the control points of the optimized curve section) can be assigned to each search space point of the fourth control point.

(110) Using this procedure, it is possible to determine, for each search space point in the discrete search spaces of the last (i.e. N-th) control point, an optimized energy value and an optimized curve (including the radius value at each control point of the optimized curve), which starts at and which ends at the respective search space point of the last control point.

(111) The data processing system then determines the search space point in the search spaces of the last control point, which maximizes the objective function OF. The corresponding curve, including the radius values for each control points of the curve is then determined to be the matched template.

(112) Thereby, a computationally efficient process is provided for determining the extent of the anatomical structure.