MULTI-MATERIAL DECOMPOSITION FOR SPECTRAL COMPUTED TOMOGRAPHY

Abstract

A computer system (MD) and relates method for spectral-data based material decomposition. The system comprises a statistical module (SM) configured to fit, per patch in input spectral imagery, a set of probability distributions (Pk) to a respective vector spectral diagram (5) for the respective patch (A). The patch is one of a plurality of patches in the input spectral imagery obtained by operation of a spectral imaging apparatus. The probability distributions are combinable into a probability map indicative of material type probabilities per image location in the input spectral imagery.

Claims

1. A computer system for spectral material decomposition, comprising: a memory that stores a plurality of instructions; and a processor that couples to the memory and is configured to execute the plurality of instructions to: fit, per patch in input spectral imagery, a set of probability distributions to a respective vector spectral diagram for the respective patch, the patch being one of a plurality of patches in the input spectral imagery, the input spectral imagery obtained by operation of a spectral imaging apparatus, the spectral imagery including at least two energy images, and combine the probability distributions into a probability map indicative of material type probabilities per image location in the input spectral imagery.

2. The system of claim 1, wherein the processor is configured to perform a combination operation to combine, per image location in the input spectral imagery, the set of probability distributions fitted to spectral diagrams for the one or more patches that include the respective image location, to obtain the probability map.

3. The system of claim 1, wherein the combination operation includes adapting contributions of the probability distributions to the probability map based on their likelihood values per image location.

4. The system of claim 1, wherein the combination operation includes averaging the probability distributions per image location.

5. The system of claim 1, wherein the processor is configured to classify the probability distributions according to material type.

6. The system of claim 5, wherein the number of different material types is at least four.

7. The system of claim 1, wherein the processor is configured to fit the probability distributions based on a maximum likelihood algorithm.

8. The system of claim 7, wherein the maximum likelihood algorithm includes the expectation-maximization algorithm.

9. The system of claim 1, wherein the number of probability distributions per spectral diagram is based on a model selection metric.

10. The system of claim 1, wherein the set of probability distributions represents a mixture model.

11. The system of claim 1, wherein the set of probability distributions comprise Gaussian distributions.

12. (canceled)

13. A computer-implemented method of spectral material decomposition, comprising: fitting, per patch in input spectral imagery, a set of probability distributions to a respective vector spectral diagram for the respective patch, the patch being one of a plurality of patches in the input spectral imagery, the input spectral imagery obtained by operation of a spectral imaging apparatus, the spectral imagery including at least two energy images; and combining the probability distributions into a probability map indicative of material type probabilities per image location in the input spectral imagery.

14. (canceled)

15. (canceled)

16. The method of claim 13, further comprising combining, per image location in the input spectral imagery, the set of probability distributions fitted to spectral diagrams for the one or more patches that include the respective image location, to obtain the probability map.

17. The method of claim 13, further comprising adapting contributions of the probability distributions to the probability map based on their likelihood values per image location.

18. The method of claim 13, further comprising averaging the probability distributions per image location.

19. The method of claim 13, further comprising classifying the probability distributions according to material type.

20. The method of claim 13, wherein the probability distributions are fitted based on a maximum likelihood algorithm.

21. The method of claim 20, wherein the maximum likelihood algorithm includes the expectation-maximization algorithm.

22. The method of claim 13, wherein the number of probability distributions per spectral diagram is based on a model selection metric.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Exemplary embodiments of the invention will now be described with reference to the following drawings, which, unless stated otherwise, are not to scale, wherein:

[0044] FIG. 1 shows a schematic block diagram of a spectral imaging arrangement;

[0045] FIG. 2 shows a block diagram of a dual-energy based material decomposition system;

[0046] FIG. 3 shows further components of the material decomposition system;

[0047] FIG. 4 shows a soft-thresholding curve as may be used in the material decomposition system;

[0048] FIG. 5 shows a flow chart of a method of spectral image-based material decomposition;

[0049] FIG. 6 shows a spectral diagram clustered based on fitted probability distributions;

[0050] FIG. 7 shows a library of probability distributions pre-classified according to material type; and

[0051] FIG. 8 shows schematic exemplary representations of probability maps obtainable by the material decomposition method or system.

DETAILED DESCRIPTION OF EMBODIMENTS

[0052] FIG. 1 schematically illustrates an example imaging arrangement IAR, including an X-ray imaging apparatus SIA and an image processing system MD. The imaging apparatus may be a CT scanner SIA. The scanner SIA includes a stationary gantry SG and a rotating gantry RG. The rotating gantry is rotatably supported by the stationary gantry SG and rotates around an examination region ER and a portion of an object or subject therein about a Z-axis. A radiation source XS, such as an X-ray tube, is supported by and rotates with the rotating gantry RG around the examination region ER. The radiation source XS emits in general wideband polychromatic X-ray radiation that is collimated to form a generally fan, wedge, or cone shaped X-ray radiation beam that traverses the examination region ER.

[0053] A radiation sensitive detector array of detector D may subtend an angular arc opposite the radiation source XS across the examination region ER. The detector array includes one or more rows of detector pixels that are arranged with respect to each other along the Z-axis and detects X-ray radiation traversing the examination region ER. The detector provides projection (raw) data.

[0054] A reconstruction system RECON, referred to herein as the reconstructor, reconstructs the projection data, generating spectral image data. As will be explored in more detail below, the system MD processes the spectral image data to obtain processed spectral image data.

[0055] A subject support PC, such as a couch, supports a subject or object (e.g., a phantom) in the examination region ER. The subject support PC is movable in coordination with performing an imaging procedure so move the subject or object with respect to the examination region ER for loading, scanning, and/or unloading the subject or object.

[0056] An operator console OC may include a human readable output device such as a display monitor DD, etc. and a user input device such as a keyboard, mouse, etc. The console OC further includes a processor (e.g., a central processing unit (CPU), a microprocessor, etc.) and computer readable storage medium (which excludes transitory medium) such as physical memory. The operator console OC allows user to control the imaging procedure.

[0057] The arrangement IAR may further include a processing unit PU such as workstation to image-process raw projection data acquired by the imager. The operator console OC and workstation may be arranged in the same computing system. The reconstructor RECON may run on the workstation PU.

[0058] Whilst the principles disclosed herein are described with main reference to CT or other volumetric/rotational imaging modalities such as C-arm imagers or similar, they are of equal application to projection imaging in radiography.

[0059] In more detail, during imaging operation, the patient, or at least a region of interest (“ROI”) of patient, resides in the examination region ER. For example, the patient may lie on the patient couch PC arranged at least partly inside the donut shaped CT examination region ER.

[0060] The X-ray source XS is energized. X-radiation emerges from the source XS, traverses the examination region ER and the ROI/patient, and is then registered at the far end at the X-ray sensitive pixels that make up an X-ray sensitive surface of the X-ray detector D.

[0061] The impinging X-radiation causes the X-ray sensitive pixels to respond with electrical signals. The electrical signals are processed by data acquisition circuitry of the scanner IA to produce digital projection raw data. The projection data in projection domain may be processed by the re-constructer RECON to compute, for example, cross sectional imagery of the ROI in image domain. The reconstructed imagery may be further processed as indeed envisaged herein into spectral imagery. The spectral imagery and/or the reconstructed imagery may be stored or displayed on a display device DD etc. The reconstructed/spectral imagery may be examined by a user for therapy or diagnostic purposes.

[0062] The re-constructer RECON is a transformer that transforms from projection domain located at the detector surface into image domain which is located in the examination region. The reconstructor RECON may implement any one or more of different types of reconstruction algorithms, such as Radon transform based algorithms, in particular filtered back-projection (FBP). Fourier-transform type algorithms are also envisaged, and so are iterative, algebraic or machine learning based reconstruction algorithms.

[0063] The reconstructed cross-sectional imagery can be thought of as image values that are assigned by the re-constructer to grid points, referred to as voxels, in the 3D portion that makes up the examination region. There may be a plurality of such sectional images in different section planes of the image domain. The plurality of section images in different such planes form a 3D image volume. Location for image values in a given cross sectional image may also be referred to herein as (image) pixels. The

[0064] As briefly mentioned, the imaging apparatus SIA, such as the CT scanner in FIG. 1, is configured for spectral imaging. The X-ray imaging apparatus SIA thus produces sets of spectral projection raw data A.sub.j, .sub.j≥.sub.2 two or more energy windows. The re-constructer RECON produces therefrom respective spectral sectional imagery I.sub.Ei.

[0065] Reconstructed spectral imagery (and the spectral projection raw data) is multi-dimensional data. This is illustrated in the upper portion in FIG. 1. In spectral imaging, for each image location L in any given image plane Π, there are associated not only one but plural image values E.sub.j, at least one such value for each energy window. For example, in dual imaging, a special instant of spectral imaging, there are two such image values E.sub.1, E.sub.2. The X-radiation emitted by the source XS is not generally monochromatic but is polychromatic. In other words, the X-ray photons that make up the X-radiation have energies in a range of values, an energy spectrum, rather than a single energy. Each such value E.sub.j for a given image location L represents contrast for different such energies or energy ranges. In spectral imaging, it is this spectral information that is harnessed. Intensities detected at the detector D depend on material and energy window, and this differentiation is carried over into image domain for the reconstructed sectional imagery. Each set j of spectral projection raw data λ.sub.j, .sub.j≥.sub.2 for each energy window E.sub.j is hence reconstructed by reconstructor RECON into a respective “energy image” I.sub.Ej, thus resulting in plural image values E.sub.j for any given image location L as mentioned above.

[0066] In this manner, spectral imaging allows resolving image contrast into plural energy windows. Resolving into two such energy windows, high E2 and low E1, is sufficient for present purposes and is usually referred to as dual energy imaging. The low energy (about < 30 keV) image I.sub.E1 represents attenuation mainly by photoelectric absorption, whilst the higher energy image I.sub.E2, (E2 > 30 keV) represents attenuation mainly by Compton scattering. Each energy image I.sub.E1, I.sub.E2 captures different aspects of attenuation.

[0067] However, resolving in more than two such energy windows, as into three, four or more such energy windows is also envisaged in spectral imaging. In embodiments, resolving into two or more than two energy values may be obtained at the detector-side by having a specially configured detector D that is equipped with counting circuitry (not shown). The counting circuitry classifies incoming intensities in projection domain into different energy bins against a set of energy thresholds. Other detector-side solutions include dual or multi-layer detector configurations. Source-side solution are also envisaged, such as dual-or multi-source imagers, or those with a single source XS equipped with fast kVp switching circuitry.

[0068] Although we will mainly refer below to dual imaging with two energy windows, E1,E2, with two energy images I.sub.E1, I.sub.E2, this is not to be construed as limiting. Specifically, principles of the present disclosure are applicable to all form of spectral imaging with more than two energy windows E.sub.j, .sub.j>.sub.2 and thus more than two energy images (for example, photon-counting based imaging) I.sub.Ej,.sub.j>.sub.2 are specifically envisaged herein in embodiments.

[0069] Spectral imaging is useful for diagnostic purposes in particular as it allows extracting more information as traditional energy integrating imaging would allow, where the spectral information described above is usually lost. In traditional imaging, the energy is integrated and is the total energy deposited that confers image contrast. The different energy images may be of interest in themselves, but the spectral imagery may be processed into other imagery still, such as a virtual monochromatic image, a contrast agent quantitative map, a virtual non-contrast image, an electron density image, and/or other spectral imagery.

[0070] Of particular interest herein is a yet further type of spectral imagery, which is referred to as material maps in material decomposition imaging. This type of imaging allows resolving image contrast into different material types. For instance, one can produce a fat image, water image and bone image and so forth, where each contrast represents the amount or concentration of each respective material fat, water, bone at a given image location L.

[0071] Applications of material decomposition imagery abound, for instance, medical image segmentation, bladder or kidney stone detection, liver-fat quantization (LFQ), and many others still.

[0072] The imaging arrangement IAR as proposed herein includes a computerized system MD configured for material decomposition, also referred to herein as the “material decompositioner MD” or simply the “decompositioner MD”.

[0073] The decompositioner MD may be implemented by one or more computing systems PU. It may be arranged centrally or on one or more servers in a distributed or cloud service architecture communicatively coupled to one or more imagers through a wired or wireless connection. Alternatively, the MD is implemented in a single computing system, such as a workstation, associated with a single imager SIA.

[0074] The spectral imagery supplied by the imager SIA may be first stored in an image storage, such as a PACS (picture archiving and communication system) in a hospital information system (HIS), and is then retrieved therefrom and processed by the decompositioner MD. Alternatively, the spectral imagery may be forwarded direct via wired or wireless communication from the imager IA or its workstation to the decompositioner MD.

[0075] The decompositioner MD may be arranged as one or more software modules that run on a single computing systems PU or, in a distributed architecture, on more than one computing systems PU, such as servers etc. A server-side arrangement MD may be useful as a group of imagers SIA in a geographical area may be supported, such as, for example, imagers SIA from different medical facilities.

[0076] The one or more computing systems PU include one or more processers, for instance capable of high-speed processing such as those of multi-core design. Parallel computing may be used. The computing system PU further includes one or more memories MEM where instructions are stored to implement the material decompositioner MD. Instead of, or in addition to, software-based implementations, at least a part of the material decompositioner MD may be arranged in hardware, for example hard-coded into dedicated circuitry, such as ASICS or others, preferably integrated into the imager SIA or its workstation.

[0077] Reference is now made to FIG. 2 which shows a schematic block diagram of the decompositioner DM including components thereof.

[0078] Input spectral imagery, including at least two different energy images, such as the high and low energy images I.sub.E1 and I.sub.E2 in dual imaging, are received at a statistical module SM.

[0079] As explained above, the input spectral imagery relates to spatial domain where each pixel/voxel represents a location (x,y,z) in the examination region ER, with respective attenuation values for different energy windows E1,E2 associated therewith. The statistical module SM as envisaged herein processes data related to the input imagery to identify plural materials present in the ROI. It is envisaged herein that the material decomposition can be done robustly for two, three, but more specifically, for four or more materials.

[0080] In more detail, the statistical model SM attempts to identify the materials based on fitting a set of probability distributions, not to the input spectral imagery themselves, but in a spectral space associated with the input spectral imagery. This spectral space, as opposed to the spatial domain in which the input spectral imagery is located, can be represented by a spectral diagram, or more specifically, by a vector spectral diagram, with dimension corresponding to the number of energy windows. For example, in dual energy imaging envisaged herein in embodiments, there are two dimensions with two axes in this spectral space indicating the respective attenuation responses in the energy windows, high and low. In photon counting applications for more than two energy windows, the spectral diagram has dimension three or higher, with three or more axes accordingly. The spectral diagram has a histogram structure, where for each point with energy coordinates (e.sub.1, e.sub.2, ... e.sub.j), there is associated the total sum of attenuation responses as per the spectral images, preferably normalized. For example, in dual imaging, the spectral diagram is thus a surface in 3D space (e.sub.1, e.sub.2) -> n, with n the normalized count over all those image locations with responses at energy combination (e.sub.1, e.sub.2).

[0081] The input energy images I.sub.E1, I.sub.E2 are hence first transformed into such a vector spectral diagram in spectral domain, and the statistical module SM operates on the said vector spectral diagram, referred in the following as the “spectral diagram”. The spectral diagram may be smoothed to facilitate the fitting operation by statistical module SM, so as to reduce the stepped appearance of the spectral diagram.

[0082] Rather than attempting to fit the probability distributions to the whole spectral diagram at once, a localized/partialized approached is envisaged herein instead as will be detailed below. The input spectral imagery is sub-divided in spatial domain into patches, sub-sets of the input spectral images. A respective spectral diagram is then prepared for each of the patches separately and the statistical module SN then fits the probability distributions per patch. A plurality of probability distributions so fitted form a respective local/partial statistical model. Because the fitting is done per patch, a set of partial models is obtained this way, one such partial model per patch. Each partial statistical model may be rendered as a partial probability map per patch.

[0083] An optional combiner Σ combines the partial models for some or all patches to obtain a global probability map Pa.sub.j for the whole image plane. In the probability map, partial or global, each entry at a given image location L represents the probability for each of the three or more materials found. The probability map may be displayed on a display device DD color- or grey-value coded.

[0084] It will be understood that the proposed method can be applied for dual or triple material decompositions, although the full benefits of the proposed scheme are most palpable in material decomposition tasks for more than three materials, such as four, five, six or even more still. Specifically, double-digit material decomposition can be done with the proposed system MD in a robust manner.

[0085] The decompositioner MD may include a use interface UI, such as a graphical user interface GUI, for the user to select the number of materials for the material decomposition to be computed. The graphical user interface (GUI) may include one or more check boxes, a dropdown menu for example. A GUI is not requirement herein, as textual input via textboxes in response to keyboard input is also envisaged herein, and so are other user interfaces UI.

[0086] Reference is now made to FIG. 3, which shows a block diagram to provide more details of the material decompositioner DM.

[0087] A patch definer PD uses one or more kernels of any shape and size to slide in a step width over image plane Π in spatial domain to define patches K1, K2 at image locations L1, L2, respectively, for spectral imagery I.sub.E1, I.sub.E2. The patches K1,2 can be of any shape and/or size and may be contiguous or may be overlapping as shown in FIG. 3 schematically. The whole image plane may be tiled with such patches K1,2, but this may not necessarily be so in embodiments where only a certain image region is analyzed. For exemplary purposes, only two patches K1,K2 at two locations L1,L2 are shown in FIG. 3, but this is for illustration only, it being understood that there are usually more than two patches depending on the size of the input imagery. The size and/or the shape of the patches can be adjusted to the data for better quality of the statistical results. On the one hand, a given patch preferably includes enough pixels to provide for better convergence of the implemented algorithm. On the other hand, it may be preferable to restrict patch size so as not to include more than the envisaged number of different types of issues and/or materials. Preferably, the size of a patch should be in the order of a few hundred pixels/voxels. In case computational resources are only of minor concern, one may use different exploratory runs of the proposed algorithm on patches of different sizes and/or shapes, and then filter out all those patches that result in partial models with only low likelihood. Another design/hyper-parameter for the proposed algorithm, apart from patch size or shape, is the amount of overlap among the patches. This amount may be quantified in terms of percentage of overlap between neighboring patches. For more robustness, a higher amount of overlapping is desirable. The design parameter for overlapping may be constant or may be adaptive. A smaller step width, and hence a higher amount of overlap, may be used for regions with strongly varying values, whilst a larger step width, and hence lower amount of overlap, may be used for more homogenous image regions.

[0088] Each such image patch K.sub.i defined by definer PD is transformed by spectral diagram generator SDG into a respective spectral diagram S.sub.i. For dual imaging, the respective local diagram S.sub.i per patch has two dimensions, one for each energy window E1,E2. Each point (e1,e2) in the spectral diagram S.sub.i in spectral space (in this case a plane) may thus represent the frequency f= n/N of measurements in the two energy images IE1,IE2 as seen at one or more spatial locations n over all locations N.

[0089] The order of operation of patch definer PD and spectral diagram generator SDG may be reversed. The whole or a part of the spectral image is transformed first into a global spectral diagram, which is then divided into a set of (partial) spectral diagrams S.sub.i. The patches in spatial domain may then be defined as the inverse image of each partial diagram S.sub.i in spatial domain. Each partial diagram S.sub.i hence gives rise to a respective patch in spatial domain, each such patch K.sub.i comprising locations whose measurement at the respective locations fall into the respective partial diagram S.sub.i. Each spectral diagram Si is hence associable with a patch K.sub.i in spatial domain.

[0090] The patches K.sub.i may or may not be topologically connected.

[0091] For each spectral diagram S.sub.i, statistical module SM builds a respective partial model by fitting parameters of a family of probability distributions to data as recorded in the respective said spectral diagram S.sub.i. A model builder MB may define a number of hyper-parameters to inform the statistical fitting operation by the statistical module SM. In embodiments, statistical mixture models are used. Such mixture models are represented by a weighted sum of individual parameterized probability densities or distributions. In embodiments, Gaussian partial mixture models GMM.sub.i are used as adapted to the respective spectral diagram S.sub.i per patch. Each probability density in the mixture model per patch Ki has its own parameters that are adapted by the statistical module, such as the respective mean and variance/standard deviation.

[0092] In particular, the number of probabilities distributions to be used for each spectral diagram S.sub.i per patch is such a hyper-parameter and this number may be defined by model builder MB in embodiments. The number of probability distributions for each partial mixture model GMM.sub.i may correspond to the number of materials one expects to find. This number may be predefined by user through user input UI or may be computed based on the spectral data by model builder MB. This number may be the same for each patch, or may vary and can be computed by model builder MB based on certain metrics as will be described in more detail above. The number of probability distributions in each partial model GMM.sub.i may be said to form a “group” or “mixture”. Whilst we refer herein to each partial mixture model as GMM.sub.i, it is understood that this partial model may not necessarily be Gaussian. Other types of probability distributions may be used instead. Also, parameterized models are not a necessity herein, and non-parametric statistical methods are also envisaged herein.

[0093] In embodiments, the statistical module SM operates to fit a linear combination of the respective group of probability distributions to the respective spectral diagram to obtain the respective partial statistical model GMM.sub.i for the given patch Ki/spectral diagram Si. The fitting operation may be done according to a statistical algorithm such as a maximum likelihood algorithm. In embodiments, to be described in more detail below, the expectation maximization algorithm is used, but other maximum likelihood or non-maximum likelihood methods are equally envisaged herein. Again, whilst Gaussian mixture models have been found to work particularly well, the present principles are not confined to Gaussian mixture modeling and other types of statistical probability distributions may be used instead. Equally, statistical mixture models, whilst beneficial herein are not necessarily the only model types that may be used herein. Non-mixture type statistical modelling approaches may be considered herein in the alternative.

[0094] The output of the statistical module SM is a respective “partial” or partial model GMM.sub.i for each patch K.sub.i/spectral diagram S¡. Because each patch K.sub.i is associated with its respective spectral diagram S.sub.i, “K.sub.i .Math. S.sub.i”, an operation on or a statement in relation to the respective (local) diagram S.sub.i may hence therefore be referred to herein as such an operation or statement per patch K.sub.i.

[0095] Operation of the statistical module may be understood as a clusterer that clusters each spectral diagram according to the number of probability distributions for the materials. However, at this stage it may still not be known which probability distribution so found represents which material. The classification of each probability distribution per patch and partial model GMM.sub.i into a spectral material type is done by a classifier CL. Each probability distribution is preferably parameterized. Parameters may include mean .Math. and variance σ.sup.2 such as is the case for Gaussian or exponential families more generally. It is these parameters that may be adjusted by the statistical module SM, based on the data in each spectral diagram per patch. The parameters describe shape and location of the respective probability distribution in the spectral diagram and hence so define the clusters. The parameters have been found to be a tell-tale of the material type. As will be explained more fully below, the classifier CL uses in embodiments a library LBR of previously stored probability distribution prototypes which are already assigned to a respective material type. The classifier CL then compares the probability distributions for similarity with the prototypes in the library to establish the material type. The parameters of the fitted and prototypical distributions may be used to determine similarity. Alternatively or instead, machine learning may be used for this classification task. The classifier CL, preferred as it is, is optional as in some cases it may only be the number of materials for each patch that are of interest for a “coarse” material decomposition analysis.

[0096] A partial model combiner Σ′ may be used to combine the partial models GMMI for a given spatial image location in the original image plane Π. It will be appreciated that in general each image location L is covered by more than one patch such is the case for overlapping patches. Specifically, each image location is thus associated not only with one, but with a plurality of such patches that cover the said location, and hence with the same plurality of spectral diagrams and, consequently, with the same plurality of respectively fitted partial models GMM.sub.i. In yet other words, each image location is thus associated with the said plurality of partial models GMMI.

[0097] It is the said plurality of partial models GMMI associated with a given image location L that are combined by the partial combiner Σ′. This partial combination operation per image location yields partial probability maps P.sub.i,j which may be of interest in themselves if only a local material decomposition analysis is sought. However, for a global analysis, a global probability map may be required. The global probability map is obtained by global combiner Σ in a further combination operation. The global combiner Σ may ensemble-average the partial probability maps P.sub.i,j to produce a global probability map Pa.sub.j as output. As will be explored in more detail below, a thresholding scheme, hard or soft, may be used in the combination operations by combiners Σ′, Σ, to down-weigh, or remove, models with low likelihood contributions.

[0098] Any one, or any sub-combination of two or more, or all of the above described components of de-compositioner MD may be integrated into a single computational/functional component. Alternatively, the decompositoner MD is implemented in distributed fashion across more than one such computational/functional components.

[0099] Operation of the above described components will now be described in more detail, using Gaussian mixture models (GMM) as an example, with the understanding that the present disclosure is applicable to any type of probability distribution with density in the exponential family, or to those with other densities.

[0100] Turning first in more detail to operation of the patch definer PD, this may be implemented by using a parametrized family of relatively small kernels of square or round shape which cover parts or the whole of the image plane, with or without overlaps. Having some or all kernels overlap is advantageous. The more overlap is applied, the more information is available for each pixel of the original image, which in turn helps secure better quality of the resulting probability map.

[0101] The family of the kernels is characterized by the size of the kernel, its shape (such as square/rectangle, curved such as circle shaped or elliptic shapes), and the step width, that is, the amount by which the kernel is shifted from a current to the next position in the image plane. It is not necessarily herein for kernel shape and/or size to remain constant as either one or both of shape and size may vary with position.

[0102] For each kernel position, a respective local spectral diagram is built by binning according to the energy density as recorded in the two (or more) energy input images IE1, IE2 for all pixel positions covered by the current kernel position.

[0103] Operation of the classifier CL based on the library of pre-classified probability distribution prototypes will now be explained in more detail.

[0104] Classification, or labelling, of the Gaussians in the respective partial models GMM.sub.i may be done by comparing the fitted Gaussians to a set of labelled pre-computed prototype Gaussians held in in library LBR, a Ground Truth Gaussians Library (GTGL), and finding one or more closet matching prototype Gaussians. This library of labelled Gaussians plays the role of the Ground Truth in the proposed algorithm. The label represents the material type (eg, fat, water, Iodine, etc) to which the respective Gaussian prototype is attributed.

[0105] The library LBR stores the set of parameters for each prototype-Gaussian, which fully defines the prototype and hence the respective material label. For example, for bivariate prototype Gaussians the set of parameters stored in the GTGL may be written as {.Math..sub.X, .Math..sub.Y, Σ.sub.XX, Σ.sub.XY, Σ.sub.YX, Σ.sub.YY}, where .Math..sub.X, .Math..sub.Y are the coordinates of the mean in the spectral plane where the respective Gaussian is centered, and Σ.sub.XX, Σ.sub.XY, Σ.sub.YX, Σ.sub.YY are the parameters of the covariance matrix that define the spread of that prototype.

[0106] The library LBR may be built up manually or in an automated fashion.

[0107] In manual labeling, the Gaussians resulting from the partial models GMMi for the respective image patches are used. Historic imagery may be used as found in PACS (picture archiving and communication system) or other image repositories. The labelling may be done by a human expert. In a good number of cases, the material label will be known a priori, such as from the nature of the historic examination. Meta-data and patient records as held alongside the imagery in the PACS etc may be used to aid labelling.

[0108] In automatic ROI labeling, a segmentation algorithm is used to segment the spatial energy images I.sub.E1, I.sub.E2 in spatial domain into known materials regions. The partial models GMM.sub.i, are then computed as described above for the corresponding spectral diagram in the spectral domain. The resulting Gaussians are then labelled accordingly and entered into the library GTGL as new prototypes. A semi-automatic version is also envisaged, where the segmentation in spatial domain is done manually.

[0109] The classification operation as administered by the classifier CL can be hard or soft. In hard classification, the single best matching prototype is found in the library LBR based on a suitable nearest neighbor measure. Interpolation in a look-up structure such as LUT may be used to interpolate over the parameters. In soft classification, the nearest neighbor match is evaluated against a threshold which may result in plural “nearest neighbor” -prototypes as candidates for the correct material label. In this soft classification embodiment, there is no single best matching probability distribution. In this case, the result is not a single material type assignment, but a set of material types with associated probabilities, each computed from the nearest neighbor candidate probability distribution from the library. The material type labels in the set are those with which the nearest neighbor candidate probability distributions from the library LBR have been labeled with.

[0110] The one or more nearest neighbors may be found by a suitable measure such as squared Euclidean distance, cross-entropy or, more generally, the Kullback-Leibler divergence:-

[00001]$\begin{array}{cc}{D}_{KL}\left(p,q\right)=ln\left(\frac{{\displaystyle .Math.{}_2}}{{\displaystyle .Math.{}_1}}\right)+\frac{{\displaystyle {.Math.}_1^2+{\left({\mu }_1-{\mu }_2\right)}^2}}{2\ast {\displaystyle .Math.{}_2^2}}-\frac{1}{2}& \text{­­­(1)}\end{array}$

wherein { .Math..sub.1, Σ.sub.1} and { .Math..sub.2, Σ.sub.2} are parameters defining the Gaussians.

[0111] As an alternative, machine leaning may be used for the classification task. A machine learning model may be trained on the above mentioned labelled ground-truth prototypes in library LBR. For example, a deep neural network with one or more hidden layer may be used, with the last layer, a softmax-layer, to provide the material type classification. The layer may be exclusive fully connected or may be fully convolutional. Alternatively, the network may include both, fully connected layers and convolutional layers. Instead of neural networks, other machine learning models may be used such as decision trees, random forests or support vector machines SVM, and other ML model configurable for classification tasks.

[0112] Turning now again to the statistical module SM, its operation based on the mixture models GMM.sub.i envisaged for the present clustering are now described in more detail. In embodiments it is assumed that the points on each spectral diagram Sj are described as a weighted mixture of n Gaussians (GMM). The density of partial model GMMi for patch K.sub.i/spectral diagram S.sub.i may be written as:

[00002]$\begin{array}{cc}f\left(x,y\right)={\displaystyle {.Math.}_{i=0}^n{p}_i}\ast N_i\left(x,y\right)& \text{­­­(2)}\end{array}$

The p.sub.i are the weights and the number “n” will be referred to herein as the “number of components”, or “component number”. The Nj are different Gaussians with different mean and/ or co-variance. The Gaussians are bi-variate for dual imaging as considered herein, but may be multi-variate for spectral imaging based on more than two energy windows. “(x,y)” is a point in the local spectral diagram Sj.

[0113] The probability of the pixel (x, y) having been generated by Gaussian G.sub.i is computed as:

[00003]$\begin{array}{cc}\text{P}\left(\left(\text{x,y}\right)\in G_i\right)=p_i\ast N_i\left(x,y\right)/{\displaystyle {.Math.}_{i=0}^n{p}_i}\ast N_i\left(x,y\right)& \text{­­­(3)}\end{array}$

The Probability of the pixel (x, y) to belong to material M.sub.j may be computed as:

[00004]$\begin{array}{cc}\text{P}\left(\left(\text{x,y}\right)\in M_j\right)=\frac{{\displaystyle .Math.{}_{N_j}{p}_j\ast N_j\left(x,y\right)}}{{\displaystyle {.Math.}_{i=0}^n{p}_i}\ast N_i\left(x,y\right)}& \text{­­­(4a)}\end{array}$

wherein N.sub.j is the set of all Gaussians attributed to material M.sub.j.

[0114] The probability distribution fitting algorithm as implemented by the statistical module SM is now described in more detail. The fitting algorithm may operate on densities of the probability distribution. The densities may be continuous or discrete. The set of the parameters defining the statistical model is optimized by improving, eg maximizing, the Likelihood or, which is mathematically equivalent, the Log-Likelihood function. The Log-Likelihood function is defined as following:

[00005]$\begin{array}{cc}L=\log p\left(X\right)={\displaystyle .Math.{}_{\left(x,y\right)}\log {\displaystyle .Math.{}_i{p}_i\ast N_i\left(x,y,{\mu }_i,{\Sigma }_i\right)}}& \text{­­­(4b)}\end{array}$

[0115] In embodiments, a Maximum-Likelihood (ML) type algorithm, such as an Expectation-Maximization (EM)-type algorithm may be used. The EM algorithm was described by A Dempster et al in “Maximum likelihood from incomplete data via the EM algorithm”, J. Roy. Stat. Soc. Ser. B 39, pp 1-38 (1977). The EM algorithm may be used for data that is generated in modes or through hidden states that generate observables. For present purposes, the different materials in the sample are the hidden states and the spectral diagrams Si are the observables.

[0116] In order to compute the set of parameters of the partial model GMM.sub.i, which maximizes the likelihood of the model for the given spectral diagram S.sub.i, an iterative process in two-steps is used, called the “Expectation” and “Maximization” step respectively. Instead of operating on the likelihood function as such, the log-likelihood function L may be used instead. Multiple iteration cycles may be run until the algorithm converges - within an agreed error margin or other stopping condition - to a maximum of the Log-likelihood function. The maximum may be a local maximum or, less likely, a global maximum.

[0117] As such, since the convergence is local, parameters computed in the standard EM setup depend on the initialization of the parameters. Also, in the standard EM algorithm, the higher the component number, the higher the values that the log-likelihood function can obtain. In the standard EM-algorithm, the component number is fixed and assumed known.

[0118] In practical modelling, the component number should be restricted. Different model selection metrics may be used for this purpose.

[0119] In embodiments envisaged herein the model builder MB computes the Bayesian information criterion (BIC) as:

[00006]$\begin{array}{cc}BIC=\ln \left(n\right)\ast k-2\ln \left(L\right)& \text{­­­(5)}\end{array}$

[0120] Alternatively, the Akaike information criterion (AIC) is used:

[00007]$\begin{array}{cc}AIC=2\ast k-2\ln \left(L\right)& \text{­­­(6)}\end{array}$

wherein: [0121] k is the number of parameters of the model GMMi; [0122] n is the number of data points in the spectral diagram for which the model is built [0123] L is the maximized value of the Likelihood/log-Likelihood of the model.

[0124] For each model GMMi for a respective patch/spectral diagram Si, a set of M Gaussians are initialized and optimized by applying the Expectation-Maximization (EM). This may be repeated in different EM-runs, for each of a number of different values for M. Specifically, for small patches K.sub.i with diameter up to about 20 pixels, the optimal value for component number is any one of M= 1, 2, 3, 4. To compute the optimal model, one computes initial models GMMi with M = {1, 2, 3, 4} by running the EM algorithm for each value of n. The final model GMMi is the one whose Likelihood/log-Likelihood function L returns the maximal value for a selection metric, such as BIC or AIC to determine the optimal model for this patch K.sub.i. BIC metric has performed well for present purpose. For each patch K.sub.i/spectral diagram S.sub.i, plural EM algorithms are run for different M’s to arrive at the respective partial model GMMi. In general, best component number M’s will depend on the patch K.sub.i.

[0125] In order to facilitate good convergence rates of each model and convergence to a good quality local maximum, it may be beneficial to select the initial centers (means) of the Gaussians at points in the spectral diagram where the density of the values are maximal or are at least higher than a threshold. Alternatively, one may place the Gaussians at points as per the means of the prototype Gaussians in the reference library LBR, if the type of materials for the decomposition are already known, as will be the case for most medical examinations/tasks. The initial covariance matrices of the Gaussians may be initialized arbitrarily or at random, since their correct values are computed in the first iterations of the EM algorithm.

[0126] By iterating the Expectation-Maximization steps, the log-likelihood L of the respective model GMM.sub.i is increasing and converges to, in general, a local maximum.

[0127] In the expectation (“E”)-step, the expectation relative to marginal densities is computed, with hidden states averaged out by summation. Probabilities of all data points (x, y) under the model GMM¡ with the current values of parameters to belong to Gaussian k is evaluated:

[00008]$\begin{array}{cc}{P}_k\left(x,y\right)=\frac{p_k\ast N_k\left(x,y,{\mu }_k,{\Sigma }_k\right)}{{\displaystyle {.Math.}_{i=0}^n{p}_i}\ast N_i\left(x,y,{\mu }_i,{\Sigma }_i\right)}& \text{­­­(7)}\end{array}$

[0128] In the maximization (M)-step, the current parameter set for all Gaussian are updated as per:-

[00009]$\begin{array}{cc}{m}_{k+1}={\displaystyle .Math.{}_{\left(x,y\right)}{P}_k\left(x,y\right)}& \text{­­­(8)}\end{array}$

[00010]$\begin{array}{cc}{p}_{k+1}=\frac{m_k}{m}& \text{­­­(9)}\end{array}$

[00011]$\begin{array}{cc}{\mu }_{k+1}=\frac{1}{m_k}{\displaystyle .Math.{}_{\left(x,y\right)}{P}_k\left(x,y\right)\ast \left(x,y\right)}& \text{­­­(10)}\end{array}$

[00012]$\begin{array}{cc}{\Sigma }_{k+1}=\frac{1}{m_k}{\displaystyle .Math.{}_{\left(x,y\right)}}{P}_k\left(x,y\right)\ast {\left(x,y-{\mu }_k\right)}^T\ast \left(x,y-{\mu }_k\right)& \text{­­­(11)}\end{array}$

In (9), “m” indicates the total number of pixels.

[0129] A motivation for the proposed local, patch-based probability distribution fitting in spectral domain, as opposed to a global fitting may be understood as follows. Many of spatially restricted fragments of the original images I.sub.E1,I.sub.E2 will represent in general one, two, or at the most three different material types. In these cases, building a partial model GMMi for a given diagram patch S.sub.i in the spectral domain will give a good statistic model, that is one with high likelihoods, which will allow good definition of probabilities for the participating materials on this patch. In particular, by using a family of overlapping patches one can obtain good enough statistics to robustly estimate material-type probabilities for all pixels of the image. The fitting operation attempts to distribute the probability mass, based on the given data in the spectral diagram. Globally, the presence of a given material may be minute, so that its contribution may be “drowned out” when attempting to fit the probability mass globally. By working locally, even small contributions can be duly accounted for. In any one sufficiently small patch, there are likely fewer materials to be found, so there is more probability mass available, and less risk for drowning out of low concentration materials.

[0130] The per-patch approach allows “fine-grained” information handling by combiners Σ, Σ′. Contributions from patches which result in models with low likelihood can be ignored in this process. Alternatively, a smooth transition from “good” fully used patches to the “bad” ignored patches can be done by using weighted averaging of the corresponding probabilities. The weighs are a function of the partial model likelihood. The global probability maps computed herein may thus be called “Ensemble Averaging Probability Maps” (EAPM).

[0131] Operation of the partial model combiner Σ′ is now described first in more detail. After the partial model GMM.sub.i has been optimized, the partial probability of a pixel (x, y) to belong to material M.sub.j may be estimated by partial model combiner Σ:

[00013]$\begin{array}{cc}\text{P}\left(\left(\text{x},\text{y}\right)\in M_j\right)=\frac{{\displaystyle .Math.{}_{N_j}{p}_j\ast N_j\left(x,y\right)}}{{\displaystyle {.Math.}_{i=0}^n{p}_i}\ast N_i\left(x,y\right)},& \text{­­­(12)}\end{array}$

wherein N.sub.j is the set of all Gaussians attributed to material M.sub.j. Eq (11) thus allows computing partial material probability map which may be sufficient for local examination.

[0132] Operation of the global combiner Σ is now described in more detail. Global combiner Σ allows combing the partial models by ensemble averaging into a global probability map.

[0133] The partial probabilities contribute to the global probabilities maps per material. We assume that for each pixel of the original image there is a plurality of local values of probabilities. We can ensure this by building the family of the kernels with center at every pixel and some non-zero radius. This approach may generate a large amount of data and hence robustness, but may increase computation time. In order to strike a decent balance, one may use a sparse family of overlapping kernels, where the distance between two centers is for example about a ⅓ of the kernel radius. Herein, sparsity relates to the amount of overlap. The sparser the less overlap there is.

[0134] In one embodiment, global combiner Σ assembles the partial probabilities into the global map by averaging per pixel. For each pixel i= (x,y) of the image, its probability to represent material M.sub.j is defined as the average of some or all probabilities defined by all partial models GMM.sub.i, whose kernel/patch K.sub.xycontains this pixel:

[00014]$\begin{array}{cc}\text{P}\left(\left(x,y\right)\in M_j\right)=\frac{{\displaystyle .Math.{{}_K}_{xy}{p}_j\left(x,y\right)}}{N\left(K_{xy}\right)}& \text{­­­(13)}\end{array}$

where N (K.sub.xy) is the number of all kernels/patches, containing this pixel (x, y). The p.sub.j’s are probabilities computed in different partial models relevant to the given pixel (x,y). Whilst for the partial maps eq (11) only the partial model for a single patch is used, in the global map eq (13) all models for patches that cover a given pixel are taken into account and averaged over, but models with likelihood may be dropped or down-weighted as will be explained in the following.

[0135] Since the mixture models envisaged herein are statistical models strongly dependent on data distribution, the quality of the models for different patches may differ widely, especially because the component number has been restricted on practical grounds. For the most of the smaller patches, up to four Gaussians will provide a reasonable Likelihood of the model, although for some patches likelihood may still deteriorate to very low values. For better results in computing the global probability map, it is proposed to apply one or two thresholds to exclude partial results with low Likelihood when combing the global probability map. With this policy, only contributions from patches with “good” statistics are used, whilst patches with “bad” statistics are ignored. This policy helps to achieve greater robustness.

[0136] For the single threshold embodiments, eq (13) above may be used. Only results from those patches, for which the partial model p.sub.j yields likelihood greater than a threshold TH, are admitted into the summation of the denominator, withp.sub.j >TH.

[0137] Rather than a hard thresholding, a weighted averaging of the partial probabilities may be used in alternative embodiments, where the weight is a function of the likelihood of the partial model:

[00015]$\begin{array}{cc}\text{P}\left(\left(\text{x},\text{y}\right)\in M_j\right)=\frac{{\displaystyle .Math.{}_{k_{xy}}C\left(L\left(k_{xy}\right)\right)p_j\left(x,y\right)}}{{\displaystyle .Math.{}_{k_{xy}}C\left(L\left(k_{xy}\right)\right)}}& \text{­­­(14)}\end{array}$

The definition of the weight function requires two thresholds T.sup.L, T.sup.H, and a transition function from zero to one. For example, a sinusoidal curve over interval (-.sub.7[, π) may be used, scaled appropriately. For example, the weight function may be defined as:

[00016]$\begin{array}{cc}C\left(L\left(k_{xy}\right)\right)=\left\{\begin{array}{c}0,ifLikelihood<LowThreshold\\ 1,ifLikelihood<HighThreshold\\ 0.5+0.5\ast \left(Sin\left(\frac{L-T^L}{T^H-T^L}-0.5\right)\ast PI\right)otherwise\end{array}\right)& \text{­­­(15)}\end{array}$

Alternatively, the soft-max-function or tanh-may be used, or others still.

[0138] FIG. 4 shows an illustration of such a smoothed threshold weight function which may be used in embodiments by the classifier CL. T.sup.L, T.sup.H indicate the low and high thresholds, respectively.

[0139] Reference is now made to the flow chart in FIG. 5 which illustrates steps of the proposed spectral-image based material decomposition method. It will be understood that the method is not necessarily tied to the architectures described above and may be understood as a teaching in its own right.

[0140] At step S510, spectral imagery is received. This may include two or more energy images. In the above described detector- or source- enabled spectral imaging modalities the energy images are natively spatially registered. In some other cases, though less preferable, a registration algorithm may be used first to spatially align the energy images.

[0141] At step S520 patches, that is sub-sets, in the spectral imagery are defined.

[0142] At step S530 for each patch in the spectral imagery, a spectral diagram is computed. Steps 520 and 530 may be reversed in order, in which case the spectral diagram is computed first for the whole spectral imagery, and the patches are defined then.

[0143] At step S540 a plurality of probability distributions are fitted per spectral diagram or patch. A mixture model such as a Gaussian mixture model may be used where each of the plural probabilities distributions are fitted in combination. A linear combination may be used. Other models may be used instead. The fitting operation results in partial models GMMi for each patch.

[0144] At an optional step S550, the fitted probability distributions in some or each partial model GMM.sub.i are classified according to material type. A library of pre-stored prototypes may be used for the classification. Machine learning may be used for classification instead.

[0145] In step S560, the partial models per image location are combined into partial probability maps, one for each image location. Combination may be done by averaging. Partial models for those patches are combined that cover that image location.

[0146] At optional step S570 the partial probability maps are averaged to obtain a global probability map for the whole image plane. The combination may include averaging, in particular ensemble averaging. The averaging may include weighted averaging of the partial models.

[0147] In embodiments, in step S570, low likelihood partial models are down-weighted relative to others with higher likelihood to so implement a soft-thresholding. A hard thresholding may be used instead, where partial models with likelihoods below a threshold are ignored entirely and are not included in the combination to form the global map.

[0148] The proposed patch-based fitting approach in spectral domain allows securing the following advantages: The local maxima may no longer depend strongly on different initializations . Even materials represented by a small number of pixels that may attract very low weights in a global model can be adequately accounted for in partial models per patch (see below the discussion of FIG. 8).

[0149] In a conventional global model fitted at once, the original position and covariance of the different Gaussians can change significantly during the optimization process, so that initial attribution of material-type may be not justifiable. But this attribution can be safely done in the proposed localized per-patch approach Gaussians obtained in a standard optimization may sometimes extend over large areas with significant overlap with other Gaussians, thus making a classification difficult, if not impossible. With per-patch fitting in spectral domain as proposed herein, such extension and overlap can be curbed. The per-patch approach allows avoiding the likelihood to grow to inflated values. If one were to apply the EM optimization to the whole image, then even the model selection metrics BIC or AIC would not help to meaningfully curb the likelihoods.

[0150] Again, whilst the above has been explained with main reference to Gaussians and Gaussian mixture models, this is not limiting and all the above is of equal application to other mixture models not necessarily involving Gaussians. Also, statistical cluster schemes other than mixture models may be used instead, such as factor analysis. Also, the groups of probability distributions {P.sup.j} of a partial model may be related other than in a linear combination as assumed in standard mixture models. More involved relationships may be used instead. Furthermore, instead of using the EM-algorithm or any of its cognates, Monte-Carlo-Markov Chain methods, moment matching, spectral methods (eg, based on singular value decomposition) or graphical methods may be used instead to fit the group of probability distributions {P.sup.j} per patch.

[0151] FIG. 6 is an illustration of a spectral diagram for the whole image, not merely for a patch. The point clouds represents points in the spectral diagram with the ellipses (four are shown) representing the clustering. Each ellipse represents a two times σ-level line of a respective Gaussian. The example model contains 4 Gaussians. Each Gaussian may be labelled as one of four materials: Air, Bone, Iodine and Kidney Stone, for example. In the global probability map, the probability of a pixel to belong to one of those materials is estimated from the fitted mixture model. Pixels on the spectral diagram may be visually encoded, for example by color-coding or grey value. The color codes for the material with highest probability. The global or partial probability maps as computed herein may be visualized on the spectral diagram (as in FIG. 6) for the whole plane of the spectral imagery, or on a patch. Alternatively, the probability map may be visualized accordingly in spatial domain, optional displayed together with the original spectral imagery, such as with one of the energy images or with a combination thereof. For example, the probability map may be superimposed as an overlay on the spectral imagery in spatial domain, with suitable color- or grey-value coding.

[0152] FIG. 7 is an illustration of a library LBR of prototype probability distributions as may be used for classification purposes as described above. Each ellipse represents a respective material type, although for each material more than one ellipse can be present in the library. The representation in FIG. 7 is in spectral domain as in FIG. 6.

[0153] FIG. 8 is an illustration of the ability of the proposed material decompositioner MD to detect even minute amounts of a material in a more than four-material type decomposition. FIG. 8 illustrates a 5-material decomposition for abdominal imagery with pane A) indicating the original spectral imagery such as the high or low image or a combination thereof. Pane B) represents a probability map for bone. Pane C) represents a probability map for soft tissue, and pane D) a probability map for Iodine. Pane E) is a probability map for fat. Notably, pane F) shows a probability map for uric acid. As can be seen, the proposed method and system MD can detect even small amounts of a material, such as in this case uric acid, indicated by the small arrow in pane F).

[0154] The proposed material decompositioner MD allows better utilizing spectral data for medical applications where automatic segmentation or material decomposition is required. The probability maps as computed herein can be used in different medical applications, such as in Liver-Fat Quantification, for obtaining virtual unenhanced images, for characterizing renal lesions, for assessing severity of liver fibrosis, to name but a few of such applications. Applications outside the medical fields are also envisaged.

[0155] Some or all of the above mentioned design/hyper-parameters may be learned by machine learning..In particular, the size and overlap percentage of the patches, the parameters of the weighting function for averaging may be determined by machine learning from a training data set.

[0156] One or more features described herein can be configured or implemented as or with circuitry encoded within a computer-readable medium, and/or combinations thereof. Circuitry may include discrete and/or integrated circuitry, a system-on-a-chip (SOC), and combinations thereof, a machine, a computer system, a processor and memory, a computer program.

[0157] In another exemplary embodiment of the present invention, a computer program or a computer program element is provided that is characterized by being adapted to execute the method steps of the method according to one of the preceding embodiments, on an appropriate system.

[0158] The computer program element might therefore be stored on a computer unit, which might also be part of an embodiment of the present invention. This computing unit may be adapted to perform or induce a performing of the steps of the method described above. Moreover, it may be adapted to operate the components of the above-described apparatus. The computing unit can be adapted to operate automatically and/or to execute the orders of a user. A computer program may be loaded into a working memory of a data processor. The data processor may thus be equipped to carry out the method of the invention.

[0159] This exemplary embodiment of the invention covers both, a computer program that right from the beginning uses the invention and a computer program that by means of an up-date turns an existing program into a program that uses the invention.

[0160] Further on, the computer program element might be able to provide all necessary steps to fulfill the procedure of an exemplary embodiment of the method as described above.

[0161] According to a further exemplary embodiment of the present invention, a computer readable medium, such as a CD-ROM, is presented wherein the computer readable medium has a computer program element stored on it which computer program element is described by the preceding section.

[0162] A computer program may be stored and/or distributed on a suitable medium (in particular, but not necessarily, a non-transitory medium), such as an optical storage medium or a solid-state medium supplied together with or as part of other hardware, but may also be distributed in other forms, such as via the internet or other wired or wireless telecommunication systems.

[0163] However, the computer program may also be presented over a network like the World Wide Web and can be downloaded into the working memory of a data processor from such a network. According to a further exemplary embodiment of the present invention, a medium for making a computer program element available for downloading is provided, which computer program element is arranged to perform a method according to one of the previously described embodiments of the invention.

[0164] It has to be noted that embodiments of the invention are described with reference to different subject matters. In particular, some embodiments are described with reference to method type claims whereas other embodiments are described with reference to the device type claims. However, a person skilled in the art will gather from the above and the following description that, unless otherwise notified, in addition to any combination of features belonging to one type of subject matter also any combination between features relating to different subject matters is considered to be disclosed with this application. However, all features can be combined providing synergetic effects that are more than the simple summation of the features.

[0165] While the invention has been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive. The invention is not limited to the disclosed embodiments. Other variations to the disclosed embodiments can be understood and effected by those skilled in the art in practicing a claimed invention, from a study of the drawings, the disclosure, and the dependent claims.

[0166] In the claims, the word “comprising” does not exclude other elements or steps, and the indefinite article “a” or “an” does not exclude a plurality. A single processor or other unit may fulfill the functions of several items re-cited in the claims. The mere fact that certain measures are re-cited in mutually different dependent claims does not indicate that a combination of these measures cannot be used to advantage. Any reference signs in the claims should not be construed as limiting the scope.