Tracking method and apparatus

Abstract

An element tracking method comprising: for each of a sequence of frames, obtaining position data comprising a respective position assigned to each of a plurality of element signal portions within said frame; and using a linking method that uses at least said assigned position data to link element signal portions represented in at least one of the frames to element signal portions represented in at least one other of the frames thereby to track movement of elements through said region of the subject, wherein at least one of a) and b): a) the linking of the element signal portions is in dependence on a measure of element density in at least part of the at least one of the frames and in at least part of the at least one other of the frames; b) the linking of the element signal portions is in dependence on a measure of element velocity in at least part of the at least one of the frames and in at least part of the at least one other of the frames.

Description

[0001] Tracking method and apparatus

[0002] Field

[0003] The present invention relates to a method and apparatus for tracking particles of contrast agent detected using medical imaging, for example, for tracking microbubbles in ultrasound imaging.

[0004] Background

[0005] Contrast-enhanced ultrasound (CELIS) is an imaging modality used in hospitals to depict the circulation of organs within the body. Microbubbles (MBs) that act as contrast agents are injected into the patient intravenously before the imaging takes place to enhance the ultrasound image. It is now available in hospitals to depict the circulation of organs within the body with improved contrast.

[0006] Super resolution ultrasound imaging (SRIII) may be obtained using particle tracking, which comprises three steps: particle detection, localisation and linking. Detection of MBs can be achieved either by using particle probability images or by using comparisons to a reference MB signal. Since they are much smaller than the imaging wavelength, they are localised, for example by identifying the centroids of the MB signals. The motion of these MBs is then tracked by linking the localised positions between consecutive frames.

[0007] Known linking methods include, for example, nearest neighbour methods and motion linking methods. Figure 1 is an illustration of a known linking methods. Figure 1 depicts a vascular structure having a first vessel and a second vessel. Figures 1 (a) and 1(b) illustrate particle tracking in accordance with two known methods. Figure 1(c) illustrates the ground truth links.

[0008] Figure 1(a) illustrates the nearest neighbour method. MBs in one CELIS frame are linked to its nearest neighbour in the next frame. For sparse MB fields, this method works well and is capable of constructing the vessel networks using MB tracks, even in low signal-to-noise ratio environments. However, as MB density increases, the number of errors produced by nearest neighbour linking also increases. To improve the linking accuracy, motion models may be added by imposing certain motion restrictions between frames, which is visualised in Figure 1 (b).

[0009] Summary

[0010] In a first aspect there is provided an element tracking method comprising: for each of a sequence of frames, obtaining position data comprising a respective position assigned to each of a plurality of element signal portions within said frame; and using a linking method that uses at least said assigned position data to link element signal portions represented in at least one of the frames to element signal portions represented in at least one other of the frames thereby to track movement of elements through said region of the subject, wherein at least one of a) and b): a) the linking of the element signal portions is in dependence on a measure of element density in at least part of the at least one of the frames and in at least part of the at least one other of the frames; b) the linking of the element signal portions is in dependence on a measure of element velocity in at least part of the at least one of the frames and in at least part of the at least one other of the frames.

[0011] The method may comprises obtaining a plurality of potential or candidate links between a first frame and a second frame of the sequence of frames and selecting one or more of the potential or candidate links in dependent on a measure of element density and / or on a measure of element velocity. The measure of element density and / or on a measure of element velocity may be of the first and / or the second frame.

[0012] The method may comprises obtaining a plurality of potential or candidate links between successive pairs of frames of the sequence of frames and, for each frame, selecting one or more of the potential or candidate links in dependence on a measure of element density and / or on a measure of element velocity. Obtaining the plurality of potential or candidate links between frames may comprise applying a motion model and / or nearest neighbour model and / or a known linking model.

[0013] The selection of candidate links may result in a plurality of links for further linking into a plurality of tracks. The selection of candidate links may be in dependence on knowledge of vascular bed properties.

[0014] The selection of one or more of the candidate links may comprise selecting links that represent or are a least indicative of a physically possible or probable movement of elements between frames. The selection of one or more of the candidate links may comprise rejecting links that represent physically impossible movement of elements between frames. A physically impossible movement may comprise movement across a physical barrier and / or between two or more physically separated channels. The selection of one or more candidate links may comprise rejecting candidate links made between different vessels, for example adjacent vessels. The linking of the element signal portions may comprise replacing a physically impossible or improbable link with a physically possible or probable link. A physically impossible or improbable link may comprise a link having a barrier and / or other restriction on its path. A physically possible link may comprise a link having no barrier or restrictions on its path. The physically impossible or improbable link may comprise a straight path and the physically possible or probable link may comprise a non-straight, for example, a curved path.

[0015] The linking of the element signal portions may comprise rejecting candidate links by applying a criteria to the candidate links based on density and / or speed and / or angle and / or velocity.

[0016] The linking of the element signal portions may result in a plurality of links, wherein each link connects a respective pair of element signal portions. The linking of the element signal portions may result in a plurality of tracks, wherein each track comprises a respective plurality of links. Each track may be representative of motion of a respective element.

[0017] Using the measure of element density and / or the measure of element velocity in the linking of the element signal portions may reduce a number of incorrect links. Using the measure of element density and / or the measure of element velocity may result in improved tracking.

[0018] The elements may also be referred to as particles. The elements may comprise contrast elements. The element signal portions may comprise contrast element signal portions.

[0019] The contrast elements may comprise microbubbles. Microbubbles are contrast enhancing agents that act as targets in ultrasound methods. A microbubble may comprise a bubble that has a size of less than one millimetre in diameter but usually larger than one micrometre. A solution of microbubbles will typically contain microbubbles that vary in size and shape. Microbubbles may have a diameter of 1 to 10 micrometres.

[0020] The element tracking method may further comprise obtaining a sequence of frames each comprising ultrasound or other medical imaging data representing an anatomical region of a human or animal subject at a respective different time, and, for each frame, identifying a plurality of element signal portions and assigning respective position data to each of the element signal portions.

[0021] Vascular knowledge driven microbubble tracking for super-resolution ultrasound imaging may be provided. An optimization of tracking microbubbles in the blood stream as they appear in ultrasound image data may be achieved. A result may be to create robust super-resolution maps of the architecture and dynamics of the vascular and microvascular bed.

[0022] The use of vascular knowledge driven microbubble tracking may contribute to accurate and early diagnosis for cancer and specifically prostate cancer patients. Currently prostate cancer is the most common cancer in men, with the second highest mortality in men and a very high unnecessary invasive intervention rate. Accurate and early diagnosis may help to address mortality and intervention rate.

[0023] A number of statistical criteria may be used that impact in the correct selection of tracks. These criteria may comprise knowledge of vascular bed properties.

[0024] The linking of the element signal portions in dependence on a measure of element density may comprise a density assisted linking method.

[0025] The method may further comprise determining the measure of element density. The measure of element density may be determined for each of a plurality of positions, for example each of a plurality of pixels. The measure of element density may be determined for each frame of the sequence of frames.

[0026] The linking of the element signal portions in dependence on a measure of element density may comprise using a density map that is representative of density in a plurality of locations, for example a plurality of pixels. The linking of the element signal portions in dependence on a measure of element density may comprise using density for all frames of the sequence of frames. The linking of the element signal portions in dependence on a measure of element density may comprise using density for a subset of the sequence of frames. The measure of density may comprise a per-pixel density.

[0027] The density assisted linking method may comprise applying at least one density criterion to potential links, for example minimum density criteria. The density assisted linking method may comprise rejecting potential links that do not meet the at least one density criterion.

[0028] The at least one density criterion may require a sufficient number of elements. The at least one density criterion may require a smooth distribution of elements.

[0029] The measure of element density may comprise an average element density. The measure of element density may comprise an average element density over pixels along a path between element signal regions. The path may be a straight line path. The at least one density criterion may comprise a minimum value for average element density along the path. A link may be rejected if an average element density along a path for that link is below the minimum value.

[0030] The measure of element density may comprise a standard deviation of element density. The measure of element density may comprise a standard deviation of element density over pixels along a path between element signal regions. The path may be a straight line path. The at least one density criterion may comprise a maximum value for standard deviation of element density along the path. A link may be rejected if a standard deviation of element density along a path for that link is above the maximum value.

[0031] The application of the continuity criteria may reject physically improbable links of an element between frames. The application of the continuity criteria may comprise defining a partial sector or other shape from the element based on speed and / or angle and determining if the link lies, at least partially, in the partial sector or other shape. The partial sector or other shape may be defined by a maximum and minimum speed and / or angle.

[0032] Use of a measure of element density may prevent, or reduce, instances of links being made between different vessels, for example adjacent vessels.

[0033] The linking of the element signal portions in dependence on a measure of element velocity may comprise a velocity assisted linking method. The velocity assisted linking method may be used to obtain a first link in a track. The velocity assisted linking method may be used to obtain at least one further link in a track.

[0034] The measure of element velocity may comprise a measure of element speed. The measure of element velocity may comprise a measure of element direction.

[0035] The method may further comprise determining the measure of element velocity. The measure of element velocity may be determined for each of a plurality of positions, for example each of a plurality of pixels. The measure of element velocity may be determined for each frame of the sequence of frames.

[0036] The velocity assisted linking method may comprise obtaining a map of speed and / or a map of direction. The map of speed and / or map of direction may be generated using previously determined tracks. The map of speed may be representative of an average speed of previously determined tracks in each of a plurality of locations, for example each of a plurality of pixels. The map of direction may be representative of an average direction of previously determined tracks in each of a plurality of locations, for example each of a plurality of pixels.

[0037] The linking of the element signal portions in dependence on the measure of element velocity may comprise using speed and / or direction for all frames of the sequence of frames. The linking of the element signal portions in dependence on a measure of element density may comprise using speed and / or direction for a subset of the sequence of frames.

[0038] The velocity assisted linking method may comprise applying at least one continuity criterion to potential links. The velocity assisted linking method may comprise rejecting potential links that do not meet the at least one continuity criterion.

[0039] The at least one continuity criterion may comprise at least one of a maximum speed, a minimum speed, a maximum angle relative to a predetermined direction, a minimum angle relative to a predetermined direction. The at least one continuity criterion may be dependent on a measure of continuity in previously determined tracks, for example a standard deviation of speed and / or a standard deviation of direction.

[0040] The linking of the element signal portions in dependence on a measure of element density may comprise a maximum density seeking method.

[0041] The maximum density seeking method may comprise determining a curved link between two element signal portions. A straight line link between the two element signal portions may previously have been rejected. The straight line link between the two element signal portions may have previously been rejected by the density assisted linking method.

[0042] The determining of the curved link may comprise finding a path with maximum element density between the two element signal portions in two consecutive frames.

[0043] The maximum density seeking method may comprise using a map of density to find a path of maximum element density. The map of density may comprise a respective element density for each of a plurality of locations, for example each of a plurality of pixels.

[0044] The maximum density seeking method may comprise adjusting a measure of velocity in accordance with a determined path, for example increasing a velocity for a path if the path is determined to be a curved path and therefore longer than a straight path. The linking method may comprise a density assisted linking method and a velocity assisted linking method. The linking method may comprise a density assisted linking method and a maximum density seeking method. The linking method may comprise a velocity assisted linking method and a maximum density seeking method. The linking method may comprise a density assisted linking method, a velocity assisted linking method and a maximum density seeking method.

[0045] The linking method may further comprise a nearest neighbour method. The nearest neighbour method may link an element signal portion in a frame of the sequence of frames to its nearest neighbour element signal portion in a next frame of the sequence. The linking method may comprise determining a set of links using a nearest neighbour method, and refining said set of links using the measure of element density and / or the measure of element velocity.

[0046] The linking method may further comprise using a motion model method. The motion model method may impose motion restrictions on linking. The linking method may comprise determining a set of links using a motion model method, and refining said set of links using the measure of element density and / or the measure of element velocity.

[0047] The linking method may comprise optimizing of a cost matrix of potential links. The optimizing of the cost matrix may comprise minimising the cost matrix.

[0048] The linking method may comprise a first run comprising a first linking method, followed by a second run comprising a second linking method. Results of the first run may be used as an input to the second run. The second linking method may be different from the first linking method. The first linking method may comprise at least one of a nearest neighbour method, a motion model method, a density assisted linking method, a velocity assisted linking method, a maximum density seeking method. The second linking method may comprise at least one of a nearest neighbour method, a motion model method, a density assisted linking method, a velocity assisted linking method, a maximum density seeking method.

[0049] The first linking method may comprise a motion model method and a density assisted linking method. The second linking method may comprise a motion model method, a velocity assisted linking method and a density assisted linking method.

[0050] The linking method may make use of neighbourhood information. Neighbourhood information may comprise information from regions of the image around and / or outside a given element signal portion, or around and / or outside of regions that are immediately adjacent to the given element. Neighbourhood information may comprise information relating to multiple elements, for example a large number of elements. Neighbourhood information may comprise information from an entire frame or multiple frames. Neighbourhood information may comprise anatomical knowledge.

[0051] At least some of the elements may be present in vessels in the human or animal subject. The method may comprise using said tracking of said movement of elements through said region to track the paths of at least some of said vessels. The method may further comprise determining a vessel map using said tracking of said movement of elements.

[0052] The linking of the element signal portions may result in a plurality of tracks, wherein each track comprises a respective plurality of links. Each track may be representative of motion of a respective element. The tracks may be representative of motion of the elements within vessels. The elements may move in accordance with blood flow. The tracks may be considered to be representative of blood flow. The tracks may be considered to provide a manifestation of vessels or of vascular structure.

[0053] At least some of the elements may be present outside vessels in the human or animal subject, and the method may comprise tracking said movement of elements outside the vessels.

[0054] The method may further comprise using said tracking to provide a measure of pulsatile motion or other motion of the subject.

[0055] The vessels may comprise blood vessels and the method may comprise using said tracking of said movement of elements through said region to track passage of blood into, out of, or through at least one anatomical feature of interest, optionally the anatomical feature comprises at least one tumour or organ.

[0056] The vessels may have a range of sizes, and the method may comprise, for at least some of the frames, identifying elements in vessels that have a range of different sizes, optionally at least some of said vessels of said range having cross-sections and / or flow rates that are at least 2 times, optionally 5 times, optionally 10 times, optionally 100 times larger than at least other of said vessels.

[0057] The determining of links and / or tracks may comprise assigning a measure of confidence to each determined link. A determining of tracks and / or vessels may be in dependence on the measure of confidence, for example omitting links and / or tracks which have a low measure of confidence. The measure of confidence may comprise a probability value.

[0058] Using a measure of confidence may result in a more robust determination of vessels. Links that are determined with low confidence may be omitted. Vessels may be determined based only on high-confidence links. Using high-confidence links may allow tracks and / or vessels to be determined with greater accuracy.

[0059] The assigning of position data to an element signal portion may comprise fitting a mathematical function to determine a position using at least one of: intensity, shape and size of signal portion. The mathematical function may comprise a Gaussian function, optionally wherein the determined position corresponds to a peak of the Gaussian function.

[0060] The method may further comprise assigning velocity or speed data to the one or more signal portions.

[0061] The identifying may comprise performing a segmentation process. The segmentation process may comprise applying at least one of a fitting, filtering and / or transform process, optionally a watershed transform process.

[0062] The identifying may comprise identifying candidate element signal portions and optionally performing a thresholding or filtering process to exclude at least some of said candidate element signal portions.

[0063] For each frame, the identifying of one or more portions of the ultrasound imaging data as being representative of an element may be at least partially unconstrained by the number of portions of ultrasound imaging data identified as being representative of an element for at least one other of the frames, optionally such that the number of elements identified can be different for different ones of frames, optionally such that the number of elements identified for at least one of the frames is substantially independent of a number of elements identified for at least one other of the frames.

[0064] The method may further comprise identifying, for at least some of the frames, each of more than 10, optionally more than 100, optionally more than 500, optionally more than 1 ,000 portions of the ultrasound imaging data per frame as being representative of respective elements. The sequence of frames may represent a measurement period having a duration in between at least one of: 1 second and 10 seconds; 10 seconds and 30 seconds; 30 seconds and 1 minute; less than 5 minutes.

[0065] The sequence of frames may comprise a frame rate in the range of 10 frames per second to 50 frames per second, optionally higher than 50 frames per second.

[0066] The sequence of frames may each comprise ultrasound imaging data obtained by ultrasound measurements on a living human or animal subject.

[0067] The linking of element signal portions thereby to track movement of elements may comprise forming track segments between consecutive frames using the position data and joining the formed track segments to produce a plurality of element tracks.

[0068] The joining of the track segments may comprise at least one of gap closing, merging and splitting, optionally based on at least one of size, distance, signal intensity and motion direction.

[0069] The method may further comprise introducing the contrast elements into the subject, optionally using at least one of bolus or continuous infusion.

[0070] The contrast elements may comprise microbubbles or any suitable contrast media, for example, nanoparticles or contrast agent particles (CAP). The other medical imaging data may comprise computed tomography (CT) data, magnetic resonance (MR) data or positron emission tomography (PET) data.

[0071] In a further aspect, which may be provided independently, there is provided an image processing system comprising a processing resource configured to: for each of a sequence of frames, obtain position data comprising a respective position assigned to each of a plurality of element signal portions within said frame; and use a linking method that uses at least said assigned position data to link contrast element signal portions represented in at least one of the frames to contrast element signal portions represented in at least one other of the frames thereby to track movement of contrast elements through said region of the subject, wherein at least one of a) and b):- a) the linking of the contrast element signal portions is in dependence on a measure of contrast element density in at least part of the at least one of the frames and in at least part of the at least one other of the frames; b) the linking of the contrast element signal portions is in dependence on a measure of contrast element velocity in at least part of the at least one of the frames and in at least part of the at least one other of the frames.

[0072] In another aspect, which may be provided independently, there is provided an imaging system comprising: an ultrasound scanner, or other scanner, configured to perform a scan of a human or animal subject to obtained a sequence of frames; and an image processing system as claimed or described herein configured to receive and process the sequence of frames to track movement of contrast elements through a region of the subject.

[0073] In a further aspect, which may be provided independently, there is provided a computer program product comprising computer-readable instructions that are executable to perform a method as claimed or described herein.

[0074] In a further aspect, which may be provided independently, there is provided a method of single particle tracking which considers the local neighbourhood dynamical and structural information of a particle when it is tracked in a medium through consecutive frames. The method may be applied to track microbubbles in contrast enhanced ultrasound.

[0075] Features of one aspect may be provided as features of any other aspects. For example, features of the system may be provided as features of the method and vice versa.

[0076] Brief Description

[0077] Various embodiments will now be described by way of example only, and with reference to the accompanying drawings, of which:

[0078] Figure 2 is a schematic diagram of an apparatus according to an embodiment;

[0079] Figure 2 is a flowchart of an element tracking method;

[0080] Figure 3 is a schematic diagram showing an element tracking method, in accordance with an embodiment;

[0081] Figure 4 is a flowchart a density assisted linking method, in accordance with an embodiment;

[0082] Figure 5 is a flowchart a velocity assisted linking method, in accordance with an embodiment;

[0083] Figure 6(a) and 6(b) illustrate the linking methods of Figure 4 and Figure 5; Figure 7 is a flowchart of a density assisted linking method, in accordance with a further embodiment;

[0084] Figure 8(a) to 8(d) is an illustration of the density assisted linking method of Figure 7;

[0085] Figure 9 is a schematic diagram of an ultrasound imaging apparatus, in accordance with an embodiment;

[0086] Figures 10(a) to 10(e) depict results obtained using a particle tracking method in accordance with an embodiment;

[0087] Figure 11 is a table of results obtained using a particle tracking method in accordance with an embodiment;

[0088] Figure 12(a) to 12(d) depict results obtained using a particle tracking method in accordance with an embodiment;

[0089] Figures 13(a) to 13(f) depict results obtained using a particle tracking method in accordance with an embodiment, and

[0090] Figure 14 is a table of results obtained using a particle tracking method in accordance with an embodiment.

[0091] Detailed Description

[0092] In general super-resolution ultrasound imaging (SRIII) can be considered to be formed of three stages: particle detection, particle localisation and linking. The following embodiments relate to methods of tracking elements comprising particles, in particular, microbubbles (MBs). Microbubbles are an example of a contrast element, Detection of MBs can be achieved using a number of known methods, for example, by using particle probability images or by using comparisons to a reference MB signal. Since the MBs are much smaller than the imaging wavelength, they need to be localised, for example, by identifying the centroids of the MB signals. Motion of the MBs is then tracked using a linking process, by linking the localised positions between consecutive frames.

[0093] Figure 2 shows a flowchart outlining the main steps of an element tracking method 10. The method 10 is directed to processing medical images of an anatomical region, for example, of a human or animal subject, following administration of a suitable contrast medium, for example microbubbles, to the subject. The medical images are captured over a period of time to allow the movement of the contrast medium through the anatomical region to be analysed. The period of time may have a duration of, for example, seconds to minutes. Signals from microbubbles can be detected using ultrasound methods.

[0094] The contrast elements may be microbubbles or any suitable contrast media, for example, nanoparticles or contrast agent particles (CAP). In other embodiments, suitable imaging techniques other than ultrasound may be used, for example, CT scans, magnetic resonance (MR) scans, positron emission tomography (PET) scans.

[0095] A first step 12 of method 10 is obtaining a plurality of frames of a frame sequence that have been acquired using ultrasound imaging. Each frame represents ultrasound image data. .

[0096] The frame sequence and image data are representative of the anatomical region, for example, of the human or animal subject, captured over a period of time. Each frame therefore comprises ultrasound imaging data that represents the anatomical region at a different time.

[0097] The sequence of frames is captured by performing ultrasound measurements on the living human or animal subject. The ultrasound measurements represent the presence of microbubbles administered to the subject. The data capture may take place at the same time as the frame processing, or data capture and frame processing may take place at different times. The sequence of frames may be stored and later obtained for the frame processing. An example ultrasound apparatus for capturing a sequence of frames is illustrated in Figure 9.

[0098] Microbubbles are contrast enhancing agents that act as targets in ultrasound methods. A microbubble is a bubble that has a typical size of less than one millimetre in diameter but larger than one micrometre. A solution of microbubbles will contain microbubbles that vary in size and shape. Typically, microbubbles have a diameter of 1 to 10 micrometres. Before and / or during the data capture stage, microbubbles are introducing into the subject using either bolus or continuous infusion.

[0099] In some embodiments, the microbubbles are infused into the subject at a rate that enables a sparse distribution of particles in the frames. A suitable infusion rate may be determined experimentally or pre-determined. Infusion rate is determined and dependent on a number of factors, for example, human physiology, microbubble suspension density, imaging system. In some embodiments, the suspension density has to be such that the microbubbles can be separated.

[0100] In some embodiments, the region to be imaged is part of the vascular bed of the subject. The microbubbles travel through the vascular bed and by imaging the microbubbles the structure of the vascular bed is revealed. A typical vessel has a size of millimetres to a few microns. In some embodiments, the vessels comprise blood vessels and the method uses tracking of said movement of microbubbles through said region to track passage of blood into, out of, or through the anatomical feature of interest. The anatomical feature may be a tumour or organ.

[0101] Each frame is composed of pixels. The pixel size is typically 100 micrometres. In some embodiments, the pixel size is in the range of about 10 to about 1000 micrometres. The pixel size is therefore larger than the microbubble to be imaged. However, the signal from the microbubble will typically have the size of several pixels in the image. This is because it will occupy the size of the point spread function (PSF). The point spread function is the response of the imaging system to the imaged microbubble.

[0102] Each frame represents a view of the region of the subject. In a given view, for example, microbubbles may appear to overlap in the obtained ultrasound image. As described elsewhere, the position resolution can be selected to be greater than the pixel resolution.

[0103] The sequence of frames may be characterised using the period of time which they represent. Any time period suitable to gather the data to image the region of interest can be used. In addition, or alternatively, the sequence of frames may be characterised using a frame rate. The frame rate may take any suitable value.

[0104] A second step 14 is directed to pre-processing the obtained sequence of frames. The second step 14 has two main components. A first component is an image registration process. This process is performed on each frame of the sequence of frames. In some embodiments, the image registration process is a rigid image registration. The image registration process acts to generate a substantially aligned sequence of frames or video loop. The registration process acts to substantially remove image deformation or image artefacts from externally induced motion, for example, operator probe movement. In other embodiments, an alternative image registration process is performed, for example, a non-rigid image registration. In some embodiments, the image registration process may be optional.

[0105] A second component of the pre-processing step, is a filtering process performed on each frame of the sequence of frames. The filtering process removes image artefacts, noise and other speckle. The filtering process may be optional, or may be performed in part, depending on the quality of the data and image.

[0106] Following pre-processing of the frames, the next step 16 is a training process. The training process allows determination of an optimized parameter set for the subsequent particle detection and classification process. The training process is performed using the same apparatus as the method and ensure than an initial optimized parameter set is used. In some embodiments, the training process is optional or may be performed separately from the other steps of the method. In some embodiments, in place of the training process or in place of part of the training process, one or more pre-determined parameters may be used. Without the training process, the remaining steps of the method may be performed to produce results. The training process may be replaced by a manual training process or manual selection of parameters. Any suitable method for selection of parameters may be used.

[0107] In some embodiments, the training process is optional. For example, the training process may not be performed and the method carried out using a pre-selected parameter set.

[0108] While shown schematically as a separate method step, in some embodiments, the training process forms part of the other steps of the method.

[0109] The training process determines an optimized parameter set. The set of parameters determined and used are described in further detail below. In order to obtain an optimized parameter set, an initial manual assessment is required. The training process may be performed on a subset of the frame sequence or on a separately provided training set of frames. The training data may be obtained separately.

[0110] At step 18, a linking process is initiated. The method 10 involves selecting a first frame from the frame sequence, and then, at step 20, identifying a plurality of element signal portions in each frame. Step 20 corresponds to identifying one or more signal portions of the ultrasound imaging data of the frame as being representative of a microbubble or plurality of microbubbles. Step 20 may also include the step of classifying the one or more identified signal portions as either being representative of a single microbubble (a single microbubble signal portion) or multiple microbubbles (a multiple microbubble signal portion).

[0111] At step 22, position data is assigned to each element signal portion. As such, position data is assigned to each of the single microbubble signal portions and to each of the multiple microbubble signal portions. In the present embodiment, the assigning of position data to an element signal portion includes fitting a mathematical function to determine a position using at least one of: intensity, shape and size of signal portion. In some embodiments, the mathematical function is a Gaussian function, and the determined position corresponds to a peak of the Gaussian function. Information about the signal portions may be stored in a memory resource, for later use by a linking model (step 22). The stored information may include, for example: particle and path position data, velocity and classification data. In some embodiments at least some of the following data is stored and used by the process: particle segmentation, particle position (e.g. localization), particle position update (e.g. with Kalman filter), particle paths (e.g. density map), particle velocity (e.g. speed map), particle motion direction (e.g. one, two or any suitable number of rose diagrams).

[0112] Following identification and assigning of position data for the first frame, the process returns to step 18, and a second frame is obtained. The identification of element signal portions and the assigning of position data is then performed on the second frame. Steps 18 and step 20 are repeated until all frames of the frame sequence, or until a pre-set number of frames of the frame sequence have undergone the detection and classification process.

[0113] Following steps 18 and 20, each frame has corresponding stored data related to a number of identified and classified single or multiple microbubble signal portions within the frame.

[0114] At the subsequent step 24, signal portions representing microbubbles across the different frames are linked together using a linking method. The linking model takes at least some of the corresponding stored data of the frame as input. The linking model uses at least said assigned position data to link single or multiple microbubble signal portions represented in at least one of the frames to single or multiple microbubble signal portions represented in at least one other of the frames. By linking signal portions across frames, the movement of microbubbles through the anatomical region of the human or animal subject is tracked. The linking of frames is described in further detail in the following.

[0115] It will be understood that steps 18 to 22 are looped over all frames in the sequence. In some embodiments, steps 18 and at least part of step 24 are looped over pairs of adjacent frames, for example, a frame and its immediate successor in the sequence, such that the identification of signal portions (step 20), assigning position data to each portion (step 22) and identification of links are looped (step 24, in part) are performed over adjacent pairs of frames in a sequence. In further detail, the sequence of frames are arranged in a time-ordered series (1 , ... N) and the steps are performed on the first pair of adjacent frames (1 and 2) to identify signal portions, assign position data and obtain links between those frames, followed by repeated those steps on the second pair of frames (2 and 3) until all adjacent pairs have been processed. Figure 3 describes an element linking method in accordance with embodiment. It will be understood that, in some embodiments, one or more steps, optionally all steps of Figure 2, are included in the method of Figure 3. For example, one or more steps of Figure 2, optionally all steps, may be performed during each run of Figure 3. In other embodiments, the method of Figure 3 is performed after alternative particle detection and localisation method. In Figure 3, a method of particle detection and subsequent localisation is represented by stage 32.

[0116] In the method of Figure 3, the MB density is measured to define criteria that each link must meet in order to be accepted. These criteria form part of a density assisted linking method, which may prevent links from being made between adjacent vessels and are described in further detail with reference to Figure 4. Secondly, the collective speed and direction information within a local neighbourhood to guide the linking. This is called velocity assisted linking, which may assist in determining any link between frames. In particular, the velocity assisted linking procedure offers advantages in circumstances where motion models cannot be used, for example, when forming the first link in a track. This is described in further detail with reference to Figure 5.

[0117] The method of Figure 3 is implemented by running the particle tracking process twice, referred to as a two run process as depicted in Figure 3. At step 33 an MB density map is obtained for use during the first and / or second run. The MB density map is described in further detail with reference to Figure 4.

[0118] In the first run (R1 , 38), MB tracking is performed using, for example, a motion model, at stage 38. As described in further detail with reference to Figure 4, the linking process includes using a step 38 of using a motion model or nearest neighbour model to obtain potential links between frames. At step 40, density assisted linking is applied to refine the tracks and reject potential links that do not satisfy density criteria. The result is a number of tracks formed by the identified links.

[0119] In the second run (R2, 36), shown in the right dashed box, a velocity assisted linking method is performed. In this embodiment, the linking process uses the velocity information gathered in run 1. This may further refine the linking after density assisted linking has been applied, particularly for the first link in a track. As described in further detail with reference to Figure 5, the velocity assisted linking method includes a step 44 of identifying potential links between frames using motion model (or other models), followed by a step 46 of refining the potential links using velocity based or continuity criteria. At step 48, a maximum density seeking method is performed to re-evaluate any links that were rejected during the density assisted linking of run 1 in order to find whether an alternate, bending link exists (at step 52). The maximum density seeking method is described in further detail with reference to Figure 7. The result of run 2 is the final track information (at step 52) including straight and curved tracks.

[0120] While Figure 3 depicts a linking method that has a first run comprising a first linking method, followed by a second run comprising a second linking method and in which results of the first run may be used as an input to the second run, it will be understood that the linking methods described in the following (the density assisted, velocity assisted and maximum density) may be performed independently or in a different order. In further embodiments, one or more linking methods may be repeated and information from a previous run may be used in a subsequent run. For example, previously obtained knowledge of the motion of blood in vessels may be used. The process may be refined over a plurality of runs using different linking methods.

[0121] Figure 4 depicts a density assisted linking method in accordance with an embodiment. As described in the following, the method links element signal portions in dependence on a measure of element density.

[0122] At steps 60 and 62, a number of candidate links are obtained using the position data of microbubbles in frames. The candidate links are formed using known methods. A non-limiting example of forming candidate links is provided in the following. Candidate links may also be referred to as potential links.

[0123] At step 60, for each identified MB in a frame a position in the next frame is predicted. As such, for a plurality of identified MBs in a first frame, plurality of predicted positions in the subsequent frame is obtained. In the present embodiment, the predicted positions are obtained using a motion model or motion continuity model. The MBs identified in frame t are labelled mt where i = 1, . . , M, and M is the number of MBs identified in frame t. The predicted positions for frame t+1 are labelled m-, again where i = 1, . , , M, and M is the number of MBs identified in frame t.

[0124] At step 62, the predicted positions for the plurality of MBs are compared to detected positions of MBs in the next frame. In the next frame (t+1) the detected positions are labelled n7- where j = and N is the number of MBs identified in frame t +1. A comparison between the predicted positions of the MBs for the next frame and the identified positions of MBs in the next frame is performed. In the present embodiments, a cost matrix is constructed and populated with measures of distance between the predicted position of each MB and the identified positions in the next frame. For M identified MBs in the first frame, and N identified MBs in the second frame, the cost matrix is an MxN matrix:

[0125] In the cost matrix, d is a measure of distance between predicted MB positions for frame t +1 and actual MB positions in frame t +1. The distance is defined as:

[0126] In the above questions, m- = (mx' , my' ) and n = (nx, ny) are the co-ordinates of the MBs. Candidate or potential links are determined by minimizing each row of the cost matrix so that each MB is linked to and by one MB in two consecutive frames.

[0127] At step 64, a density map is obtained. The density map is representative of density in a plurality of locations, for example a plurality of pixels. In the present embodiment, the density map represents a value of a measure of density, in this embodiment, an average density, over all frame for each pixel. To obtain the map, a value of average MB density is calculated for each pixel and stored in a density map, also referred to as a pixel density map. The average MB density is the number of MBs localised to this pixel for all image frames. The subsequent density criteria may use the density map or alternative measures of density over all or a subset of frames. In some embodiments, the density map may be obtained at another stage of the method or determined prior to step 60 of Figure 4.

[0128] At step 66, density assisted linking criteria (also referred to as simply density criteria) are applied to the candidate links. In the present embodiment, the criteria are applied to all possible links represented in the cost matrix. These criteria prevent links from being made between different vessels. To apply the density criteria, a measure of microbubble (MB) density is determined for a plurality of positions. In the present embodiment, the measure of MB density is calculated for each of a plurality of pixels and for each frame in the sequence of frames. Density assisted linking may prevent links from being made between different vessels. For an image sequence under investigation, the MB density (x,y), in pixel (x,y) is defined as the number of MBs localised in this pixel for all image frames. Since MBs travel within vessels, correct links between subsequent frames (n^ and n;) should satisfy minimum MB density criteria, which in turn update the cost matrix. This update of the cost matrix ensures that links are rejected if they do not satisfy the density criteria.

[0129] The first density criteria uses an average element density over pixels along a path between pixels or other signal regions. In the present embodiment, an average MB density is calculated along a path of a candidate link. Such a measure is calculated by averaging over value of MB density along the path. In the present embodiment, these paths are straight line paths. The first density criteria is applied by calculating an average MB density for each candidate link and compared the average MB density to a minimum value. If the average MB density is below the minimum value the link is rejected at step 68. In the present embodiment, rejecting the link comprises updating the cost matrix by assigning an infinite value to the cost matrix.

[0130] In the present embodiment, the first density criteria is represented as follows: ’

[0131] 0) where p is the average MB density over all the pixels along the straight-line path between mtand n;. Equation 3 ensures that there is sufficient MB density along the trajectory between the start and end pixels of a link, which will prevent links from being made between neighbouring vessels. The parameter cpis a control parameter for accepting or rejecting a link. In this embodiment, the control parameter is defined as a percentage + The control parameter may be determined from a previous run. In some embodiments, the control parameter may be refined over a number of runs. In some embodiments, the control parameter is dependent on a statistical measure obtained from previous runs, for example, a standard deviation or a multiple of a standard deviation, for example, 1.5 of the standard deviation.

[0132] A second density criteria used a measure of MB density in the form of standard deviation of a MB density over all pixels along a straight line path. In the present embodiment, the second density criteria is applied by calculating the average standard deviation for each candidate link and comparing the comparing the average standard deviation to a maximum value for standard deviation. If the average standard deviation is above the maximum value, the link is rejected at step 68. In the present embodiment, rejecting the link comprises updating the cost matrix by assigning an infinite value to the cost matrix.

[0133] In the present embodiment, the second density criteria is represented as follows: where o-pis the standard deviation of the MB density over all the pixels along the straight-line path between mtand n;and cais a control parameter. Equation 4 ensures that the MB density does not fluctuate significantly along the path between mt and n;. The parameter cacontrols how large the fluctuation of the MB density can be before the link is rejected. The parameter camay be selected based on previous runs. For example, a value for the parameter may be selected (for example, the value may be 1 / 3 or other suitable fraction) and during a previous run and the track density variations obtained during that run may be used to determine standard deviations of all tracks. This information may be then be used to select the parameter. The control parameter may be determined from a previous run. In some embodiments, the control parameter may be refined over a number of runs. In some embodiments, the control parameter is set as the standard deviation or a multiple of the standard deviation.

[0134] The first criteria restricts links between MB to satisfy a minimum MB density criteria. The second criteria ensures that the MB density does not fluctuate significantly along a track. Use of a measure of MB density may prevent, or reduce, instances of links being made between different vessels, for example adjacent vessels. The density criteria ensure that along the trajectory of a track, there must be a sufficient number and a smooth distribution of MBs, thereby avoiding jumps of MBs between vessels and improving the vessel reconstruction.

[0135] At step 69, a number of tracks are determined using the identified links. Each track will be understood as being formed from a number of links. Once formed, the track is representative of motion of a respective element through vessels, for example, in accordance with blood flow.

[0136] Figure 5 depicts a velocity assisted linking method in accordance with an embodiment. As described in the following, the method links element signal portions in dependence on a measure of velocity. At steps 70 and 72 a motion model is used to identify potential links between frames. Step 70 and 72 correspond to steps 60 and 62 of Figure 4.

[0137] At step 74, speed and direction maps are obtained. In the present embodiment, the speed and direction maps are obtained using previously determined track information. For example, as depicted in Figure 3, the speed and direction maps are obtained during a second run. That second run follows a first run in which track information is obtained and that track information is used to determine the speed and direction maps. All of the determined MB tracks are used to generate the speed and direction maps. In other embodiments, the tracks may be obtained using alternative methods otherthan density assisted linking. In general, while Figure 3 depicts a first and second, either one or both of the runs may be iteratively performed using information obtained from a previous run. As such, the obtained tracks may be refined over a number of iterations.

[0138] The map of speed is representative of an average speed of previously determined tracks in each of a plurality of locations, for example each of a plurality of pixels. In the present embodiment, an average speed of all tracks is calculated for each pixel. The pixel values of average speed form the speed map. In addition to a map of average speed, a map of the standard deviation of speed of all tracks is also generated at this stage. The map of standard deviation of speed is determined by calculating the standard deviation of speed of all tracks passing through each pixel o-r(x,y).

[0139] At step 74, a direction map is also determined using track information. The direction map is representative of an average direction of previously determined tracks in each of a plurality of locations, for example, each of a plurality of pixels. In the present embodiment, an average direction for all tracks is calculated for each pixel. The pixel values of average direction form the direction map. In addition to a map of average direction, a map of the standard deviation of direction of all tracks is also generated at this stage. The map of standard deviation of direction is determined by calculating the standard deviation of direction of all tracks passing through each pixel ag(x,y).

[0140] At step 76, velocity based criteria are applied to the candidate links identified at step 72. The velocity based criteria may also be referred to as continuity criteria. When a MB travels along a vessel, its speed and direction vary between neighbouring frames but should be within a certain continuity range. In this embodiment, two velocity assistance criteria, fixed and variable continuities respectively are applied. The criteria that provides the best results may be dependent on the results obtained in R1. The fixed continuity criterion may provide better results in the circumstances in which all the tracks in the entire image have fairly good continuity in both speed and direction, and their variations both locally and globally are small and well defined. In some embodiments, the application of the continuity criteria may represent rejection of physically improbably or impossible paths of a microbubble.

[0141] In the present embodiment, the first criteria on link speed and direction is applied. The application of the speed criteria includes first comparing the link speed to a lower bound (corresponding to a minimum speed) and to an upper bound (corresponding to a maximum speed). The application of the direction criteria includes comparing the link direction to a maximum angle relative to a predetermined direction and a minimum angle relative to a predetermine direction. If the speed of the link does is not in the range defined by the maximum and minimum speeds or the direction is not in the range defined by the maximum and minimum angles, the link is rejected.

[0142] The first criteria may be represented as follows: where v and 6 are the speed and direction of the link between MBs mt and n;, = (mx, my) are the pixel coordinates of MB nt; in frame t. The coefficients and p2are the scaling parameters giving the fixed lower and upper boundaries of the speed respectively, where / ?i < 1 and p2> 1. The angle 0 defines the fixed angle variation allowed around 6 (mi) In the present embodiment, all three parameters are estimated from the first run. For example, after the first run, an angle difference and maximum angle deviations can be calculated from existing tracks to give an approximate value. In some embodiments, the criteria may be represented as a shape, for example, a partial annulus or a partial sector, or an arc having an thickness defined by the upper and lower bounds on speed. In an embodiment, the angle subtended by the partial sector is 0. The first, inner radius of the partial sector is the lower bound on speed in Equation (5) and the second, outer radius of the partial sector is the upper bound on speed in Equation (5). The application of the criteria is determining if the link lies inside the defined shape. In some embodiments, the parameters of the criteria may be determined from previous runs and / or are refined over a number of runs.

[0143] In the present embodiment, a second criteria on link speed and direction is applied. The second criteria applies when the speed and direction continuity is not consistent. In the present embodiment, a variability criteria is applied based on a measure of continuity in previously determined tracks. In the present embodiment, a standard deviation of speed and direction are determined for each track, for each pixel, and used to applied the criteria. The second criteria includes comparing a measure of link velocity to upper and lower bounds determined from per-pixel velocity and corresponding standard deviation per-pixel.

[0144] The second criteria is represented as:

[0145] In the above criteria, y is a parameter that controls the narrowness of the continuity range. A lower value results in a more restrictive continuity range. As seen in Equation (6), the restriction varies from pixel to pixel depending on the standard deviation of the speed and direction in that pixel. The second criterion may be applied when the speed and direction continuity is not consistent, varying significantly across the image. In an embodiment, if variations are not uniform across an image, the method may use local variations. In some embodiments, statistics for a previous run are measured at more than one pixel. If a determination is made that the variations are sufficiently uniform across an image based on a previous run, then a maximum and minimum of speed and angle may be set for all pixels. On the contrary, if a determination is made that the variations are not uniform across an image based on a previous run, then these may be set at a pixel level.

[0146] As described above, links that are rejecting by failing to satisfy either criterion. In the present embodiment, rejecting the link comprises updating the cost matrix by assigning an infinite value to the cost matrix. At step 80, tracks are formed using the cost matrix, for example, as described with reference to step 69 of Figure 4.

[0147] Figure 6(a) illustrates the concept of density assisted linking and Figure 6(b) illustrates the concept of velocity assisted linking. Figure 6(a) and 6(b) both illustrate a vascular structure 102 having a first vessel 104 and a second vessel 106. The grey arrows depict direction of blood flow through the first and second vessels. Figure 6(a) and 6(b) illustrates microbubbles in different frames.

[0148] Figure 6(a) illustrates a MB 108 in a first frame (at time, t) and four MBs in a second, subsequent frame (at time t+1). The four MBs at the second subsequent frame are referred to as: first MB 110a, second MB 110b, third MB 110c and fourth MB 110d. Four candidate (or potential) links are depicted in Figure 6(a) between each respective MB in the second frame and the MB of the first frame. These are determined, for example, using a motion based or nearest neighbour model, as described above. Four such potential links (112a, 112b, 112c, 112d) are depicted in Figure 6(a).

[0149] Likewise, Figure 6(b) illustrates a first MB 114 in a first frame (at time, t) and four MBs in a second, subsequent frame (at time t+1). The four MBs in the second frame are labelled: 116a, 116b, 116c and 116d. The velocity assisted linking method restricts the number of potential links to a single potential link 118 between the MB 114 in the first frame and the third MB 116c of the second frame. A shape in the form of a partial sector 119 is drawn from the first microbubble. The partial sector is an arc or partial annulus. The sector has a subtended angle defined by the allowable angle range above. The thickness of the sector (between the inner and outer radius) is the defined by the lower and upper bounds on the speed. It can be seen that the shape allows physically probable links (for example, link 118) to be selected from candidate links.

[0150] Figure 7 is a flowchart of a maximum density method in accordance with an embodiment. As described with reference to Figure 4, a linking method, for example, a density assisted linking method may reject straight line links between positions. If a straight line link is accepted then it can be used to form part of a track. If the straight line link is rejected it may be wrong. Alternatively, the link may have the correct start and end points but may not be straight and may have a non-straight trajectory. A maximum density based linking method is described in the following. It will be understood that the method of Figure 7 is performed with reference to a pixel map, as described above. At step 82, the start and end positions for a rejected link are obtained. In the present embodiment these correspond to an initial and final pixel, labelled mtand n;. At step 84, for the initial pixel, a neighbourhood window is applied about the initial pixel. The neighbourhood is a 3x3 array of pixels thereby defining 8 neighbours.

[0151] At step 86, one of the eight neighbouring pixels is selected based on at least their MB density values, for example, as obtained from the pixel map. In the present embodiment, the neighbouring pixels are filtered based on a threshold density value, such that only pixels above the threshold are considered. In addition, the filtered pixel that is closest in distance to the end point of the link is selected as the next pixel on the new path.

[0152] At step 88, a new neighbourhood is defined about the pixel selected at step 86 and a further pixel selected. At this stage, the window is slid to a new position and centred about the selected pixel to define eight new neighbours. The selection of one of the neighbouring pixels is performed as in step 88 based by filtering pixels and selecting a filtered pixel closest in distance to the end point. The further selected pixel is added as the next pixel on the new path.

[0153] At step 90, the selection of pixels on the new path is repeated until the final pixel is reached. Once the final pixel is reach, a new path is defined through the series of selected pixels. The new path has a non-straight trajectory. Each identified pixel is labelled in a sequence of (%1;y ... (xN,yN). If the method does not result in (xN,yN) reaching the end pixel the method terminates and the link is rejected as no alternative path is possible.

[0154] Figure 8 illustrates the steps of the maximum density method. Figure 8(a) depicts a pixel map in which the density of each pixel value is represented by a shade. For the purpose of Figure 8(a) to 8(d), the pixel densities are indicated on a scale of nine possible values from zero density (first density 202) to highest density (ninth density 218): first density 202, second density 204, third density 206, fourth density 208, fifth density 210, sixth density 212, seventh density 214, eight density 216, ninth density 218. While Figure 8 depicts densities on a scale, it will be understood that these labels are provided for illustration only, and the density may run on a continuous scale. The same density labels of Figure 8(a) apply to Figure 8(b) to (d) but are omitted for clarity. The pixels can be divided into higher density pixels (corresponding to labels 210, 212, 214, 216 and 218), lower density pixels (correspond to labels 204, 206, 208) and zero density pixels (corresponding to label 202). Zero density means that substantially no MBs are detected. The black circle 220 represents the starting pixel and the black triangle 222 represents the ending pixel. The black circle corresponds to mtand the black triangle corresponds to n;. Figure 8(a) depicts a first candidate link: the straight-line link 224. As part of the maximum density method, this candidate link is rejected because it violates the density criteria. The straight-line link is from mtand n;and is rejected by the path density criteria

[0155] Figure 8(b) depicts inspection in the neighbourhood around the starting pixel 220 (corresponding to step 84). The neighbourhood is a 3x3 pixel area 230 centred on the starting pixel 222. Of the eight neighbouring pixels in pixel area 230, the pixel (above a pre-determined threshold) that is nearest to the end pixel 220 is identified as the first identified pixel 228 of the new path and is assigned co-ordinates (%i,yi).

[0156] At Figure 8(c) the process is repeated by inspection in a further neighbourhood around the first identified pixel 228, corresponding to step 86. A further neighbourhood is defined a further 3x3 pixel area 232 centred on the first identified pixel 228. The nearest pixel to the final pixel 220 that is above the pixel threshold is identified as second identified pixel 234 and is assigned co-ordinates (x2,y2).

[0157] This process is repeated until, as depicted in Figure 8(d) the following pixels are identified: first identified pixel 228, second identified pixel 234, third identified pixel 236, fourth identified pixel 238, fifth identified pixel 240 and sixth identified pixel 242. These are labelled as (x^y- ... (%6<ye)- In this example, the method terminates because (x6,y6) is equal to the final pixel 222. A curved or bending link 244 is then fitted between the starting pixel 220 to the final pixel 222 via the identified pixels to form the new link.

[0158] The determining of the curved link between a pre-determined initial and final part can be regarded as finding a path based on density values and proximity to the end of the rejected. Such a path may comprise adjusting a measure of velocity in accordance with a determined path, for example increasing a velocity for a path if the path is determined to be a curved path and therefore longer than a straight path. The straight link may represent a path that is physically impossible, in that a restriction or other barrier is on the path. The curved link replacing the straight link may represent a path that is physically possible, in that the path is free of restrictions or barriers or is a path. For example, the rejected straight link may represent a path between different vessels or a path over a vessel boundary while the curved link may represent a path contained inside a single vessel. An embodiment of an ultrasound imaging apparatus 100 is illustrated schematically in Figure 9. A linear array transducer 110 is configured to transmit energy into an object to be imaged (for example, a part of the human or animal body) and to receive ultrasound echoes from the object. In other embodiments, any suitable transducer may be used to receive 2D or 3D ultrasound data and may not be a linear array. For example, the transducer may be a phased array, curvilinear array, or other suitable array.

[0159] The received ultrasound echoes are digitized by an analogue to digital converter (ADC) 112. The digitized ultrasound data is stored in memory 116.

[0160] The digitized ultrasound data is processed by a processor 114 and the resulting processed data may also be stored in memory 116, or in another memory. The processor 114 may be configured to perform beamforming of the digitized echo data. In some embodiments, more than one processor 114 may be used.

[0161] The ultrasound imaging apparatus 100 also includes a screen 118 for the display of ultrasound images (which results from further post-processing that is tailored to the display requirements) and one or more user input devices 120 (for example, a keyboard, mouse or trackball) for receiving input from a user of the ultrasound imaging apparatus 112, for example a sonographer.

[0162] In the embodiment of Figure 9 the ultrasound imaging apparatus 100 is configured to obtain ultrasound data using the linear array transducer 110 and to process that data. In other embodiments, a separate processing apparatus (for example, a workstation or general purpose computer) may be used to process ultrasound data that has previously been acquired by an ultrasound machine.

[0163] In some embodiments, the organs imaged by the may include any human or animal organ, for example: organs of the digestive system including, for example, liver, pancreas; organs of the urinary system including, for example, kidneys or bladder; organs of the cardiovascular system, including the heart; any sensory organ; organ of the central nervous system, including the brain. The anatomical feature may include any component of the lymphatic system, including, for example, lymphatic vessels and lymph nodes. The anatomical feature may include any component of the cardiovascular system, including, for example, the heart, arteries, veins, capillaries, a part of the vascular bed. As the microbubbles move inside a vessel, a first application, as described above, is to track the microbubbles to delineate the vessels. However, there is also vessel movement due to the movement of the subject, for example, due to heart beating and breathing. The tracks therefore contain information, not only on blood movement but also on vessel wall movement due to the pulsatile motion of the vessel as propagated by the heart. A further application is therefore to use the tracks to perform a measurement or to represent pulsatile motion of the subject.

[0164] In some embodiments, for example, in CT, MR or PET applications, in contrast to microbubbles, the contrast agents may move outside the vascular space. Therefore another application is to track movement of the contrast agents outside the vessels. This movement may be slower.

[0165] In the above-described embodiments, a step of identifying signal portions is described. Further comments on a non-limiting example of identifying microbubble or element signal portions are provided in the following. As a first example, a binary image from an original greyscale image is generated. This step involves generating multi-scale Haar-like features which measure local contrast in different shapes and sizes. Haar-like features are formed using non-local and statistic mapping of the original greyscale image and leads to the binary image. Using the generated binary image, one or more signal portions are identified in the generated binary image that are representative of a microbubble or plurality of microbubbles. Signal portions that represent microbubbles of every size and intensity are identified. The signal portions identified at this stage are candidate signal portions. Some candidate signal portions may correspond to noise. There is no constraint imposed on the number of signal portions that are identified. Therefore, all candidate signal regions are identified at this stage irrespective of their properties. A final number of signal regions, and therefore the number of microbubbles represented, is determined by later filtering and selection stages. A first number of signal portions / microbubbles identified in a first frame can be different to a second number of signal portions / microbubbles identified in a second frame, and the first number and the second number can be independent.

[0166] A first filtering and / or thresholding process is then applied to the generated binary image that acts to exclude signal portions that do not correspond to microbubbles or do not meet other criteria. The first filtering process filters any detected pixels that are isolated. A strict or less strict connectivity filter is applied. The first filtering process also includes a thresholding process. The size of signal region is determined and compared to a pre-determined threshold value. If the size of the signal region is less that a minimum size, the signal region is discarded. The intensity of the signal region is also determined and compared to a pre-determined threshold value. The intensity of the signal region may be determined, for example, by summing the intensity of all the pixels of the signal region. If the intensity of the signal region is lower than the threshold intensity value, the signal region is discarded. The filtering and / or thresholding process may not be performed, or may be performed only in part, depending on the quality of the data and image.

[0167] A particle probability image (PPI) is then generated from the binary image. Each detected signal portion is enhanced. Each signal portion of the binary image can be refined based on foreground and background values of the PPI. An image smoothing process is then applied to the original grayscale image. The smoothing process includes using a Gaussian smoothing kernel on the original image to provide a convolved image. In addition, local maxima are found and refined using the following criteria: the local maxima must be in a segmented region and have a PPI value must be above a pre-determined particle region threshold. A second filtering process is applied. The second filtering process includes generating a watershed transform and discarding any detected signal portions that are too large. A position is then assigned to each remaining signal portion. A geometric weighted centroid is determined for each signal portion using both the size and intensity of the signal portion. Values of size, shape and intensity may be determined at this step or may be re-used from a previous step, for example, those determined at step 36. Position may be assigned at a resolution that is higher than the resolution of the original image. An optional step, is classifying the identified signal portion as corresponding to a single or multiple microbubble event. The classification may be based on the size, intensity and shape of the signal portion. Whilst signal portions corresponding to multiple microbubbles may be assigned only a single position, they may remain classified as multiple microbubble events to enable multiple paths in the later linking stage and assign a particle density that is closer to the correct one.

[0168] In addition, in the above-describe embodiments, a motion linking model is described. The following description of a motion model are provided, however, it will be understood that alternative linking using other motion models or nearest neighbour models may be used. In alternate embodiments, a different motion model is implemented as part of the energy matrix updating. For example, a non-linear motion model may be used based on the assumption that the microbubbles represented by the signal portions move according to a non-linear motion model.

[0169] A first stage of the linking step using a motion model involves linking of identified single or multiple microbubble signal portions to track movement of microbubbles between consecutive frames to form track segments using the position data for the single and multiple microbubble signal portions. The second stage involves joining the formed track segments to produce a plurality of microbubble tracks. The first stage of linking signal portions frame to frame to form track segments is based on a parameter, the maximum displacement, labelled MD, which controls the number of pixels which are permitted for a particle to move from a frame to a subsequent frame in the sequence. The first stage includes calculating an energy matrix based on a cost analysis process. The energy matrix is determined using a nearest neighbour method between signal portions of consecutive frames. Each neighbourhood is defined by a value MD, which is the maximum allowable displacement between frames. The energy matrix may include information such as particle intensity, size and shape in order to facilitate the recognition of a particle in its next location. The linking model is operable to link at least one single microbubble signal portion in at least one of the frames to a multiple microbubble signal portion in a subsequent at least one other of the frames or vice versa. The linking model is also operable to link multiple microbubble signal portion representing a different number of microbubbles in different frames. Following the initial linking stage, the energy matrix is optimized and updated using a multiple Kalman filter and the assumption that the microbubbles represented by the signal portions move according to a linear motion model.

[0170] The following non-limiting comments on experimental results are provided. To test the new linking strategies described above (also referred to as PTNS - Particle Tracking with Neighbourhood Similarities), the methods are applied to three different data sets: synthetic, animal and human prostate CELIS images. These data sets will allow us to quantitatively assess the performance of PTNS under various different conditions, some with distinct and known characteristics that can be compared to the structural and dynamical features recovered by PTNS.

[0171] A synthetic flow model is first generated to simulate blood flow in a vascular network. Pointlike particles are then injected into the network, moving along the flow. These particles are then blurred using a variable point spread function that mimics MBs observed in real in vivo CELIS data. White Gaussian noise is finally added to each frame. This leads to a realistic CELIS data set of MBs travelling with variable speeds along a pre-designed network structure, each of which is tagged with an ID number to be used as the ground truth to compare with the tracking results of PTNS.

[0172] To evaluate the performance of PTNS, the speed map of MBs obtained by different particle tracking methods are compared to the ground truth, which is shown in Figure 10. The speed map registers the average speed of all the tracks that pass through each pixel in the image sequence which, as shown, displays not only the dynamics of the flow but also the structure of the vessel network; here the network has a symmetric structure, with high speeds at the top and bottom of the network and low speeds in the middle. The ground truth is shown in Figure 10(a).

[0173] The nearest neighbour linking method combined with the motion model is applied to the synthetic data. As shown in Figure 10(b) the overall structure of the network is mostly present but some of the structure details are missing. While the speed in the middle of the network is close to the ground truth, it is significantly lower at both ends. This is because the nearest neighbour method tends to link nearby MBs, and as such this method, even with motion models, can lead to significant linking errors, particularly in high speed environments.

[0174] PTNS is applied to the synthetic data. The results for R1 are shown in Figure 10(c), which has shown a significant improvement in the recovery of both the network structure and the flow dynamics by comparing with the ground truth. This must be attributable to the density assisted linking as it is the main difference to the nearest neighbour method. However, there are still some noticeable wrong links in R1 , mainly in the middle region of Figure 10(c) where tracks are more densely populated. Since the results of R1 show that the velocity continuity is mostly uniform across the network, velocity assisted linking model equation 5 can be applied to perform the second run.

[0175] As shown in Figure 10(d), most of the wrong links in the middle region in R1 are removed, due to the linking being based on the velocity information gathered in R1. The maximum density seeking method is applied to include bending links, the results of which is shown in Figure 10(e). While there does not appear visually to be much difference compared to 10(d), the bending links do improve the results, which can be seen through the quantitative analysis in Table 1.

[0176] Figure 11 depicts a table (T able 1 ) that shows the number of correct links that a linking method produced (CL), the total number of links produced by the method (DL) and the total number of links in the ground truth (GT). From these, two statistics are calculated: precision and Jaccard index, to quantitatively measure the performance of different linking methods. From left to right in the table: the nearest neighbour plus motion model method, PTNS R1 , PTNS R2 with straight-line linking, PTNS R2 with bending linking, PTNS R2 with combined straight-line and bending linking. It is observed that the precision and Jaccard index improves as PTNS is applied. On applying the maximum density seeking method in R2, there is a slight decrease in precision but a large increase in the Jaccard index. PTNS therefore not only improves the tracking accuracy but also uses the data more efficiently, both of which are important in a clinical setting.

[0177] PTNS is also tested on an animal data set. The data was acquired from a sheep ovary using a 1.2mL bolus injection of SonoVue (Bracco, Geneva, Switzerland) contrast agent. 936 CELIS frames are collected with a frame rate of 5 Hz. Given that this is an in vivo data set, there is no ground truth. However, a good judgement about the performance of PTNS can be made due to the simplicity of the data. There are also very striking features which will test the ability of PTNS to reconstruct structures using MB tracks. As such, PTNS is applied to this data set as described with reference to Figure 3. The results are shown in Figure 12.

[0178] The MBs in each CELIS frame have their positions localised and collected into a MB number map, where each pixel (x, y) in the MB number contains the total number of MBs localised in that pixel for all image frames, as previously described in section 2.2. This result is shown in Figure 12(a). Within the MB number map, there are three distinct features that can clearly be seen. These are labelled 1 , 2 and 3 in Figure 12(a) and they are as follows: a horseshoe shaped structure, a curved structure, and a thin vessel. Reproduction of these three features will be used in place of the ground truth for this data set. Just as with the synthetic data, speed maps are used to show the dynamics and structure, which allows a clear visualisation of the performance of each tracking method.

[0179] The nearest neighbour linking method is applied to the sheep data. From Figure 12(b) it can clearly be seen that feature 1 is not reproduced well. The horseshoe shape is heavily polluted, which shows that there are many links being made across the gap in the middle of the structure. Feature 2 is mostly well reproduced, but there are links connecting it to feature 1 on the left hand side which are incorrect. Of the three features, feature 3 is the best reproduced. In this region, the MB concentration is low, the speed is low and the structure is thin, which are all conditions that the nearest neighbour method performs well in.

[0180] Next PTNS is applied to sheep data. The speed map for R2 using only straight-line links is shown in Figure 12(c), and the speed map for R2 combining straight-line and bending links is shown in Figure 12(d). Features 1 and 2 are better reproduced using PTNS than they are using the nearest neighbour method. The empty region in the middle of the horseshoe shape in feature 1 is preserved in Figure 12(c), and there are no longer any incorrect links connecting the left of feature 2 to feature 1. However, straight-line links are not enough to preserve the connection between features 1 and 2 on the right side of feature 2. Combining straight-line links with bending links reproduces this connection. It also further improves the structure of features 1 and 2, which can be seen in Figure 12(d).

[0181] The thin vessel in feature 3 is not well reproduced in Figure 12(c). There is a slight improvement in Figure 12(d), though there are still parts missing. Because the MB density is sparse in this region, it is more likely for links to be rejected due to density assisted linking. Because of this sparsity of MBs, there isn't a clear path of highest MB density for the maximum density seeking method to follow, and thus it fails to find alternative, non-straight-line paths. As a result, this shows that the maximum density seeking method is not suitable for data sets with a very low MB concentration.

[0182] PTNS is then applied to human prostate data. CELIS data was collected in the Western General Hospital in Edinburgh from patients scheduled for radical prostatectomy with full ethical approval. An infusion of contrast agent was administered to the patient with a roughly constant flow rate, and data was collected for 3-4 minutes at a frame rate of 10 Hz using an iU22 Philips scanner with C10-3v transducer. As with the sheep data set, there is no ground truth, but certain structural and dynamical features and post-pathology results can be used to indicate the performance of the tracking method. The results for this test are shown in Figure 13.

[0183] Figure 13(a) is the B-mode prostate image of a patient under investigation, where the prostate border is circled with a yellow line. The MB number map is shown in figure 13(b) by collecting all the MBs from the 2070 CELIS frames of the prostate, the MB number density patterns in the map indicating the vascular network, though without dynamical information. PTNS is applied to track these MBs, beginning with R1. The track number map for R1 , which is the number of tracks passing through each pixel, is shown in figure 13(c) while the corresponding speed map is shown in figure 13(d). The cancer region is identified by the pathology team and is encircled by the red dashed line in figure 13(c). As seen, the cancer region has both a high track number and a high speed. It is noted that a high track number and high speed is also found in healthy areas, noticeably in the central region. However, there are no clear spatial structures in the cancer region compared to those in healthy areas.

[0184] R2 of PTNS is then applied, with a variable continuity range as defined by equation 6. Results with different continuity levels have been investigated by adjusting the control parameter y. Figures 13(e) and 13(f) are the track number map and the speed map for R2 with y = 1.5. Comparing figure 13(e) with figure 13(c), it can clearly be seen that the total number of tracks in the prostate is reduced in R2, which also shows clearer structures due to the imposed continuity restrictions. Similar features due to continuity restrictions are also observed when comparing figure 13(f) with figure 13(d).

[0185] From the literature, it is known that the vasculature of prostate cancer is different than that of healthy tissue (Forster et al., 2017; Alizadeh et al., 2013; Vau-pel and Kelleher, 2012). Cancer blood vessels are typically thicker than healthy blood vessels (Forster et al., 2017; Muller et al., 2008), with increasingly tortuous vessels (Alizadeh et al., 2013) and a heterogeneous blood flow (Vaupel and Kelleher, 2012; Jochumsen et al., 2020). As such, some discrepancies in the structures and dynamics may be expected between the cancer region and healthy region when the continuity range is varied. In order to quantify this, a comparison between the number of tracks in the cancer and healthy regions at different levels of continuity restrictiveness is investigated.

[0186] Figure 14 includes a table (Table 2) compares the track number in the cancer and healthy regions of the patient (SRI010) for R1 , and R2 with two different y-vlaues. As seen in table 2, while both the cancer and healthy region shows a decrease in the number of tracks as the continuity range becomes more restrictive, the cancer region shows a greater reduction in track number; 5% for y = 1 .5 and 9% for y = 1.

[0187] As such, the cancer region shows a larger reduction in track number than the healthy region, meaning that the former lacks motion continuity more than the latter. This in turn indicates that the vessels and corresponding blood flow in cancer regions behave more chaotically in structure and dynamics respectively. This finding is in line with literature observations of prostate cancer blood vessels (Alizadeh et al., 2013; Vaupel and Kelleher, 2012; Jochumsen et al., 2020). These findings show that PTNS can potentially be a useful tool to identify distinctive structural and dynamical features in cancer areas to assist in the diagnosis of prostate cancer.

[0188] A new method of Particle Tracking with Neighbourhood Similarities which seeks to address some important issues in current particle linking approaches and methods when applied to CELIS images. In synthetic data, it has been shown that PTNS improves the performance of particle tracking. By increasing the precision and Jacqard index, PTNS uses the data more efficiently than previous methods. In animal data, it has been demonstrated that PTNS can construct complex structures in in vivo data despite lacking ground truth. Finally, PTNS is applied to human prostate data to investigate the structure and dynamics of the vasculature for a patient with prostate cancer. Because PTNS is a two run process, the second run may be used to probe different levels of motion continuity by varying the control parameters, which can provide the means to distinguish differences in behaviour between cancer and healthy tissue. Using this approach, it is observed that the cancer region exhibited a more chaotic vascular structure, resulting in a larger decrease in track number when the velocity continuity was made more restrictive. This is consistent with the current understanding of prostate cancer. More broadly, this result shows that super resolution ultrasound imaging of prostate cancer may potentially be developed into tools for diagnosis and focal therapy.

[0189] Although description of particular embodiments has been provided above, it should be understood that these embodiments are illustrative only and that the claims are not limited to those embodiments. Those skilled in the art will be able to make modifications and alternatives to the described embodiments which are contemplated as falling within the scope of the appended claims. Each feature disclosed or illustrated in the present specification may be incorporated in any embodiment, whether alone or in any appropriate combination with any other feature disclosed or illustrated herein. In particular, one of ordinary skill in the art will understand that one or more of the features of the embodiments of the present disclosure described above with reference to the drawings may produce effects or provide advantages when used in isolation from one or more of the other features of the embodiments of the present disclosure and that different combinations of the features are possible other than the specific combinations of the features of the embodiments of the present disclosure described above.

Claims

CLAIMS:1 . An element tracking method comprising: for each of a sequence of frames, obtaining position data comprising a respective position assigned to each of a plurality of element signal portions within said frame; and using a linking method that uses at least said assigned position data to link element signal portions represented in at least one of the frames to element signal portions represented in at least one other of the frames thereby to track movement of elements through said region of the subject, wherein at least one of a) and b):- a) the linking of the element signal portions is in dependence on a measure of element density in at least part of the at least one of the frames and in at least part of the at least one other of the frames; b) the linking of the element signal portions is in dependence on a measure of element velocity in at least part of the at least one of the frames and in at least part of the at least one other of the frames.

2. The method of claim 1 wherein the linking of the element signal portions results in a plurality of links, wherein each link connects a respective pair of element signal portions.

3. The method of any preceding claim, wherein the method may comprises obtaining a plurality of potential or candidate links between a first frame and a second frame of the sequence of frames and selecting one or more of the potential or candidate links in dependent on a measure of element density and / or on a measure of element velocity of the first and / or second frame.

4. The method of any preceding claim, wherein obtaining the plurality of potential or candidate links comprises applying a motion model and / or nearest neighbour model and / or a known linking model and / or wherein the selection of one or more candidate links comprises applying a criteria to the potential or candidate links to obtain physically possible and / or probable links between frames.

5. The method of any preceding claim, wherein the linking of the element signal portions may result in a plurality of tracks, wherein each track comprises a respective plurality of links, optionally wherein each track is representative of motion of a respective element.

6. The method of any preceding claim, wherein the elements comprise particles and / or contrast elements, for example, microbubbles, optionally, wherein the microbubbles comprise a bubble between 1 micrometre in diameter to one micrometre, optionally between 1 micrometre and 10 micrometres.

7. The method any preceding claim comprising obtaining a sequence of frames each comprising ultrasound or other medical imaging data representing an anatomical region of a human or animal subject at a respective different time, and, for each frame, identifying a plurality of element signal portions and assigning respective position data to each of the element signal portions.

8. The method of any preceding claim wherein the linking of the element signal portions is in dependence on a measure of element density and wherein the method comprises determining the measure of element density, optionally, for each of a plurality of positions, for example each of a plurality of pixels and / or for each frame of the sequence of frames.

9. The method of claim 7 or 8, wherein the linking of the element signal portions in dependence on a measure of element density comprises at least one of: a) using a density map that is representative of density in a plurality of locations, for example a plurality of pixels; b) using density for all frames of the sequence of frames; c) using density for a subset of the sequence of frames;10. The method of any preceding claim, wherein the measure of element density comprises at least one of: a per-pixel density an average element density; an average element density over pixels along a path between element signal regions; a standard deviation of element density and / or a standard deviation of element density over pixels along a path between element signal regions.11 . The method of any preceding claim, wherein the method comprises a density assisted linking method comprising: applying at least one density criterion to potential links, for example minimum density criteria, optionally rejecting potential links that do not meet the at least one density criterion.

12. The method of claim 11 , wherein the at least one density criterion comprises minimum value for average element density along the path, for example, a straight line path, optionally wherein a link is rejected if an average element density along a path for that link is below the minimum value13. The method of claims 11 or 12, wherein the at least one density criterion comprises a maximum value for standard deviation of element density along the path, for example a straight line path, optionally, wherein a link may be rejected if a standard deviation of element density along a path for that link is above the maximum value.

14. The method of any preceding claim, wherein the measure of element velocity comprises a measure of element speed and / or a measure of element direction.

15. The method of any preceding claim, further comprise determining the measure of element velocity, optionally wherein the measure of element velocity is determined for each of a plurality of positions, for example each of a plurality of pixels and / or for each frame of the sequence of frames.

16. The method of any preceding claim, wherein the linking of the element signal portions in dependence on a measure of element velocity comprises a velocity assisted linking method.

17. The method of claim 16, wherein the velocity assisted linking method comprises obtaining a map of speed and / or a map of direction.

18. The method of claim 17, wherein at least one of: a) the map of speed and / or map of direction is generated using previously determined tracks; b) the map of speed is representative of an average speed of previously determined tracks in each of a plurality of locations, for example each of a plurality of pixels.; c) the map of direction is representative of an average direction of previously determined tracks in each of a plurality of locations, for example each of a plurality of pixels.

19. The method of claims 16 to 18, wherein the velocity assisted linking method comprises applying at least one continuity criterion to potential links, optionally rejecting potential links that do not meet the at least one continuity criterion.

20. The method of claim 19, wherein the at least one continuity criterion comprises at least one of: a maximum speed, a minimum speed, a maximum angle relative to a predetermined direction, a minimum angle relative to a predetermined direction and / or is dependent on a measure of continuity in previously determined tracks, for example a standard deviation of speed and / or a standard deviation of direction.21 . The method of claims 19 to 20, wherein applying the at least one continuity criterion to potential links comprises defining a partial sector or other shape from the element based on at least speed and / or angle and determining if the link lies, at least partially, in the partial sector or other shape.

22. The method of any preceding claim, wherein the linking of the element signal portions in dependence on the measure of element velocity comprise using speed and / or direction for a subset and / or for all frames of the sequence of frames.

23. The method of any preceding claim, wherein the linking of the element signal portions in dependence on a measure of element density may comprise a maximum density seeking method, optionally wherein the maximum density seeking method may comprise determining a curved link between two element signal portions, further optionally wherein the determining of the curved link comprises finding a path with maximum element density between the two element signal portions in two consecutive frames.

24. The method of claims 23, wherein the maximum density seeking method comprises using a map of density to find a path of maximum element density, optionally, wherein the map of density may comprise a respective element density for each of a plurality of locations, for example each of a plurality of pixels.

25. The method of claims 23 to 24, wherein the maximum density seeking method comprises adjusting a measure of velocity in accordance with a determined path, for example increasing a velocity for a path if the path is determined to be a curved path and therefore longer than a straight path.

26. The method of any preceding claim, wherein the linking of the element signal portions comprises replacing a physically impossible or improbable link, for example, a straight link, with a physically possible or probable link, for example, a non-straight link.

27. The method of any preceding claim, wherein the linking method further comprises a nearest neighbour method and / or a motion model method, wherein the method comprisesdetermining a set of links using a nearest neighbour method and / or motion model method and refining said set of links using the measure of element density and / or the measure of element velocity.

28. The method of any preceding claim, wherein at least some of the elements are present in vessels in the human or animal subject and wherein the method comprises using said tracking of said movement of elements through said region to track the paths of at least some of said vessels.

29. The method of any preceding claim, wherein the linking of the element signal portions results in a plurality of tracks, wherein each track comprises a respective plurality of links.

30. An image processing system comprising a processing resource configured to: for each of a sequence of frames, obtain position data comprising a respective position assigned to each of a plurality of element signal portions within said frame; and use a linking method that uses at least said assigned position data to link contrast element signal portions represented in at least one of the frames to contrast element signal portions represented in at least one other of the frames thereby to track movement of contrast elements through said region of the subject, wherein at least one of a) and b):- a) the linking of the contrast element signal portions is in dependence on a measure of contrast element density in at least part of the at least one of the frames and in at least part of the at least one other of the frames; b) the linking of the contrast element signal portions is in dependence on a measure of contrast element velocity in at least part of the at least one of the frames and in at least part of the at least one other of the frames.

31. An imaging system comprising: an ultrasound scanner, or other scanner, configured to perform a scan of a human or animal subject to obtained a sequence of frames; and an image processing system as claimed or described herein configured to receive and process the sequence of frames to track movement of contrast elements through a region of the subject.

32. A computer program product comprising computer-readable instructions that are executable to perform the method of claims 1 to 29.

Images