Tracking method and apparatus

Abstract

An element tracking method comprising the steps of: obtaining, for each successive frame, position data including a respective position assigned to each of a plurality of element signal portions within the frame; and a linking method using at least the 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 tracking element movement through the region of the subject, the linking method using at least one of a) and b): a) the linking of the element signal portions depends on a measure of element density in at least a portion of at least one of the frames and at least one other at least portion of the frames; and b) the linking of the element signal portions depends on a measure of element velocity in at least a portion of at least one of the frames and at least one other at least portion of the frames.

Description

[Technical Field]

[0001] The present invention relates to a method and apparatus for tracking contrast agent particles detected using medical images, for example for tracking microbubbles in ultrasound images. [Background technology]

[0002] Contrast-enhanced ultrasound (CEUS) is an imaging modality used in hospitals to visualize the circulation of internal organs. A contrast agent, called microbubbles (MB), is injected intravenously into the patient before imaging, enhancing the ultrasound image. It is currently used in hospitals to visualize the circulation of internal organs with enhanced contrast.

[0003] Super-resolution ultrasound imaging (SRUI) can be obtained using particle tracking techniques. MB detection can be achieved using particle probability images or comparison with reference MB signals. Because MBs are much smaller than the imaging wavelength, they can be localized by identifying the center of the MB signal. MB movement can then be tracked by linking their positional relationships between successive frames. Summary of the Invention [Problem to be solved by the invention]

[0004] Known linking methods include, for example, nearest neighbor and motion linking methods. Figure 1 illustrates known linking methods. Figure 1 shows a vasculature having a first vessel and a second vessel. Figures 1(a) and 1(b) illustrate particle tracking by two known methods. Figure 1(c) illustrates ground truth linking.

[0005] Figure 1(a) shows the nearest neighbor method. MBs in one CEUS frame are linked to the nearest MBs in the next frame. When MBs are sparse, this method works well and can construct a vascular network using MB trajectories even in low signal-to-noise environments. However, as MB density increases, the number of errors introduced by nearest-neighbor links also increases. To improve link accuracy, a motion model can be added by imposing certain motion constraints between frames. [Means for solving the problem]

[0006] In a first aspect, there is provided a method for tracking an element, comprising the steps of: acquiring, for each successive frame, position data including a respective position assigned to each of a plurality of element signal portions in a frame; and using the acquired position data to link the element signal portions in the one frame to element signal portions in other frames in order to track movement of the element through an area of ​​a subject, the element tracking method comprising at least one of a) and b) below: a) the linking of the elementary signal portions depends on a measure of elementary density in at least a portion of the one frame and in at least a portion of the other frame; or or b) the linking of said elementary signal portions depends on measurements of elementary velocities in at least part of said one frame and in at least part of said other frame;

[0007] The method may further include obtaining a plurality of potential or candidate links between a first frame and a second frame of the consecutive frames, and selecting one or more of the potential or candidate links depending on element density measurements and / or element velocity measurements of the first and / or second frames.

[0008] The method may include obtaining a plurality of potential or candidate links between pairs of consecutive frames, and for each frame, selecting one or more of the potential or candidate links depending on a measure of element density and / or a measure of element velocity. Obtaining a plurality of potential or candidate links between frames may include applying a motion model and / or a nearest neighbor model and / or a known link model.

[0009] The selection of candidate links may form multiple links for further linking to multiple tracks. The selection of candidate links may depend on knowledge of the characteristics of the vascular bed.

[0010] The selection of one or more of the candidate links may include selecting links that represent, or at least indicate, physically possible or probable movement of the element between the frames. The selection of one or more of the candidate links may include rejecting links that represent physically impossible movement of the element between the frames. Physically impossible movement may include movement across a physical barrier and / or movement between two or more physically separated channels. The selection of two or more candidate links may include rejecting candidate links made between different vessels, for example, adjacent vessels.

[0011] The linking of the component signal portions may include replacing physically impossible or impossible links with physically possible or plausible links. Physically impossible or impossible links may include links that have barriers and / or other restrictions on their paths. Physically possible links consist of links that do not have barriers or restrictions on their paths. Physically impossible or impossible links consist of straight paths, and physically possible or plausible links consist of non-straight, e.g., curved, paths.

[0012] The linking of the elementary signal portions may include rejecting candidate links by applying density and / or velocity and / or angle and / or speed based criteria to the candidate links.

[0013] The element signal portion links may form a plurality of links, each link linking a respective pair of element signal portions. The element signal portion links may form a plurality of tracks, each track including a respective plurality of links. Each track may represent the motion of a respective element.

[0014] The use of element density measurements and / or element velocity measurements in the element signal portion links can reduce the number of erroneous links. The use of element density measurements and / or element velocity measurements can improve tracking.

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

[0016] The imaging element may include microbubbles, which are contrast-enhancing agents that act as targets in ultrasound imaging. Microbubbles may include gas bubbles with diameters less than 1 millimeter, but typically larger than 1 micrometer. Microbubble solutions typically contain microbubbles of various sizes and shapes. Microbubbles have diameters greater than 1 micrometer and smaller than 10 micrometers.

[0017] The element tracking method may further include the steps of acquiring successive frames each containing ultrasound or other medical image data representing an anatomical region of a human or animal subject at a 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.

[0018] Vascular knowledge-driven microbubble tracking for super-resolution ultrasound imaging may be provided. Optimization for tracking microbubbles in the bloodstream appearing in ultrasound image data may be achieved. As a result, robust super-resolution maps of vascular and microvascular bed structure and dynamics may be generated.

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

[0020] A number of statistical criteria are used to correctly select the tracks, including knowledge of the characteristics of the vascular bed.

[0021] Linking of element signal portions in dependence on a measure of element density can constitute a density-assisted linking method.

[0022] The method may further include determining a measure of element density. A measure of element density may be determined for each of a plurality of locations, e.g., each of a plurality of pixels. A measure of element density may be determined for each frame of successive frames.

[0023] Linking element signal portions depending on a measure of element density may include using a density map representing density at multiple locations, e.g., multiple pixels. Linking element signal portions depending on a measure of element density may include using density for all frames of the consecutive frames. Linking element signal portions depending on a measure of element density may include using density for a subset of the consecutive frames. The density measure may include density per pixel.

[0024] The density-assisted linking method may include applying at least one density criterion to potential links, such as a minimum density criterion. The density-assisted linking method may include rejecting link candidates that do not satisfy the at least one density criterion.

[0025] At least one density criterion may require a sufficient number of elements. At least one density criterion may require a smooth distribution of the elements.

[0026] The element density measure may include an average element density. The element density measure may include an average element density over pixels along a path between element signal regions. The path may be a straight path. The at least one density criterion may include a minimum value of the average element density along the path. If the average element density along the link's path is below the minimum value, the link may be rejected.

[0027] The element density measure may include a standard deviation of the element density. The element density measure may include a standard deviation of the element density across pixels along a path between element signal regions. The path may be a straight path. The at least one density criterion may include a maximum value of the standard deviation of the element density along the path. A link may be rejected if the standard deviation of the element density along the link's path is greater than or equal to the maximum value.

[0028] Applying the continuity criterion may reject physically impossible links between elements between frames. Applying the continuity criterion may include defining a partial sector or other shape from the elements based on velocity and / or angle, and determining whether the link is at least partially within the partial sector or other shape. The partial sector or other shape may be defined by a maximum and minimum velocity and / or angle.

[0029] The element density measure can be used to prevent or reduce instances of links being made between different vessels, for example, adjacent vessels.

[0030] The linking of element signal portions depending on the element velocity indicator may constitute a velocity-assisted linking method, which may be used to obtain the first link of the trajectory, and which may be used to obtain at least one further link in the trajectory.

[0031] The element velocity measurements may include element velocity measurements. The element velocity measurements may include element orientation measurements.

[0032] The method may further comprise determining a measure of element velocity. The measure of element velocity may be determined for each of a plurality of locations, e.g., for each of a plurality of pixels. The measure of element velocity may be determined for each frame of successive frames.

[0033] The velocity-assisted link method may include obtaining a velocity map and / or a direction map. The velocity map and / or direction map may be generated using previously determined tracks. The velocity map may represent an average velocity of the previously determined tracks at each of a plurality of locations, e.g., at each of a plurality of pixels. The direction map may represent an average direction of the previously determined tracks at each of a plurality of locations, e.g., at each of a plurality of pixels.

[0034] Linking element signal portions in dependence on a measure of element velocity may include using velocity and / or direction for all frames of the consecutive frames. Linking element signal portions in dependence on a measure of element density may include using velocity and / or direction for a subset of the consecutive frames.

[0035] The speed-assisted linking method may include applying at least one continuity criterion to potential links. The speed-assisted linking method may include rejecting link candidates that do not satisfy the at least one continuity criterion.

[0036] The at least one continuity criterion may include 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, etc. The at least one continuity criterion may depend on a previously determined measure of continuity in the track, such as a standard deviation of speed and / or a standard deviation of direction.

[0037] The step of linking component signal portions depending on a measure of component density may constitute a maximum density pursuit method.

[0038] The method for determining maximum density may include determining a curvilinear link between two component signal portions. A straight line link between two component signal portions may have previously been rejected. A straight line link between two component signal portions may have previously been rejected by the density-assisted linking method.

[0039] Determining the curvilinear link may involve finding a path that maximizes the element density between two element signal portions in two consecutive frames.

[0040] The method for determining the maximum density may include using a density map to find a path of maximum element density. The density map may include a respective element density for each of a plurality of locations, e.g., for each of a plurality of pixels.

[0041] The maximum density search method may include adjusting a measure of speed along the determined path, for example, increasing the speed along the path if the path is determined to be curved and therefore longer than a straight path.

[0042] The linking method may include a density-assisted linking method and a velocity-assisted linking method. The linking method may include a density-assisted linking method and a maximum density seek method. The linking method may include a velocity-assisted linking method and a maximum density seek method. The linking method may include a density-assisted linking method, a velocity-assisted linking method, and a maximum density seek method.

[0043] The linking method may further include a nearest neighbor method, which may link element signal portions in a frame of consecutive frames to their nearest element signal portions in the next frame in the sequence. The linking method may include determining a set of links using a nearest neighbor method and narrowing down the set of links using element density measurements and / or element velocity measurements.

[0044] The linking method may further include using a motion modeling method, which may impose motion constraints on the links. The linking method may include determining a set of links using the motion modeling method, and narrowing down the set of links using a measure of element density and / or a measure of element velocity.

[0045] The linking method may include optimizing a cost matrix of potential links, wherein optimizing the cost matrix may include minimizing the cost matrix.

[0046] The linking method may include a first execution including the first linking method followed by a second execution including the second linking method. Results of the first execution can be used as input for the second execution. The second linking method may be different from the first linking method. The first linking method may comprise at least one of a nearest neighbor method, a motion model method, a density-assisted linking method, a velocity-assisted linking method, and a maximum density seek method. The second linking method comprises at least one of a nearest neighbor method, a motion model method, a density-assisted linking method, a velocity-assisted linking method, and a maximum density seek method.

[0047] The first linking method may include a motion model method and a density-assisted linking method, and the second linking method may include a motion model method, a velocity-assisted linking method, and a density-assisted linking method.

[0048] The linking method may utilize neighborhood information. The neighborhood information may consist of information from regions of the image around and / or outside a given element signal portion, or around and / or outside a region immediately adjacent to a given element. The neighborhood information may include information relating to multiple elements, e.g., multiple elements. The neighborhood information may consist of information from an entire frame or multiple frames. The neighborhood information may include anatomical knowledge.

[0049] At least some of the elements may be present within a vessel within a human or animal subject. The method may include using the tracking of the movement of the elements through the region to track a path of at least a portion of the vessel. The method may further include using the tracking of the movement of the elements to determine a vascular map.

[0050] The links of the element signal portions may form a plurality of tracks, each track including a respective plurality of links. Each track may represent the movement of a respective element. The tracks may represent the movement of the element within a vessel. The elements move according to the blood flow. The tracks may represent the blood flow. The trajectory may provide a representation of the vessel or vasculature.

[0051] At least some of the elements may be present extravascularly in the human or animal subject. The method may include tracking said movement of the elements extravascularly.

[0052] The method may further include using the tracking to provide a measure of pulsatile or other movement of the subject.

[0053] The vasculature may include a blood vessel. The method may include using said tracking of movement of said element through said region to track blood flow in, out, or through at least one anatomical feature of interest (optionally, the anatomical feature includes at least one tumor or organ).

[0054] The vessels may have a range of sizes, and the method may include identifying, for at least some of the frames, elements within the vessels having a range of sizes, optionally at least some of the vessels in said range having a cross section and / or flow rate that is at least 2 times, optionally 5 times, optionally 10 times, optionally 100 times larger than at least other of the vessels.

[0055] The step of determining the links and / or tracks may include the step of assigning a reliability measure to each determined link. The determination of the tracks and / or vessels may depend on the reliability measure, for example, links and / or tracks with a low reliability measure may be omitted. The reliability measure may include a probability value.

[0056] Using a confidence measure, the determination of vessels may be more robust. Links that are determined to be less confident may be omitted. Vessels may be determined based only on highly confident links. Using highly confident links allows for greater accuracy in determining tracks and / or vessels.

[0057] Assigning the position data to the component signal portions may include fitting a mathematical function to determine the position using at least one of the intensity, shape, and size of the signal portions, wherein the mathematical function comprises a Gaussian function, and optionally the determined position corresponds to a peak of the Gaussian function.

[0058] The method may further include assigning velocity or velocity data to one or more element signal portions.

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

[0060] The identifying may include identifying candidate element signal portions and, optionally, performing a thresholding or filtering process to exclude at least some of the candidate element signal portions.

[0061] The step of identifying, for each frame, one or more portions of the ultrasound imaging data as representing elements may not be at least partially constrained by the number of portions of the ultrasound imaging data identified as representing elements for at least one other of the frames, and optionally the number of elements identified for at least one of the frames may be substantially independent of the number of elements identified for at least one other of the frames, such that the number of elements identified may be different for different ones of the frames.

[0062] The method may further include identifying, for at least some of the frames, each of more than 10 portions, optionally more than 100 portions, optionally more than 500 portions, optionally more than 1000 portions of the ultrasound imaging data per frame as representing a respective element.

[0063] Consecutive frames can represent a measurement period having a duration of at least one between, for example, 1 second and 10 seconds, 10 seconds and 30 seconds, 30 seconds and 1 minute, less than 5 minutes, etc.

[0064] The successive frames may include a frame rate in the range of 10 frames per second to 50 frames per second, optionally greater than 50 frames per second.

[0065] Each successive frame may comprise ultrasound image data obtained by ultrasound measurement on a living human or animal subject.

[0066] The step of linking element signal portions to track element movement may thereby include the steps of forming track segments between successive frames using the position data, and linking the formed track segments to generate a plurality of element tracks.

[0067] Combining of track segments may include at least one of gap closing, joining and splitting, and may optionally be based on at least one of size, distance, signal strength and direction of movement.

[0068] The method may further optionally include introducing an imaging element into the subject using at least one of a bolus injection or a continuous infusion.

[0069] The imaging elements may include microbubbles or any suitable contrast agent, such as nanoparticles or contrast agent particles (CAPs). Other medical image data may include computed tomography (CT) data, magnetic resonance (MR) data, or positron emission tomography (PET) data.

[0070] In a further aspect, which may be provided independently, there is provided an image processing system including a processing resource configured as follows. This system is obtaining, for each of the successive frames, position data including a respective position assigned to each of a plurality of element signal portions within a frame; using a linking method to link element signal portions in one frame to element signal portions in other frames using the acquired position data to track element movement through an area of ​​the subject; The processing resource includes a processing resource configured to: This system has the following characteristics a) or b). a) The linking of the elementary signal portions depends on measurements of elementary density in at least part of the one frame and in at least part of the other frame. b) the linking of said elementary signal portions depends on measurements of elementary velocities in at least part of said one frame and in at least part of said other frame;

[0071] In another aspect, which may be provided independently, 1. An imaging system comprising: an ultrasound or other scanner configured to perform a scan of a human or animal subject to acquire successive frames; an image processing system as previously described configured to receive and process the successive frames to track movement of the imaging element through a region of the subject; an imaging system comprising: is provided.

[0072] In a further aspect, which may be provided independently, there is provided a computer program product comprising computer readable instructions operable to carry out the methods described herein.

[0073] In a further aspect, which may be provided independently, a method for single particle tracking is provided that takes into account the mechanical and structural information of the local neighborhood of a particle as it is tracked in a medium through successive frames, which method can be applied to tracking microbubbles in contrast-enhanced ultrasound.

[0074] Features of one aspect may be provided as features of another aspect, e.g., system features may be provided as method features, and vice versa. [Brief explanation of the drawings]

[0075] Various embodiments will now be described, by way of example only, with reference to the accompanying drawings, in which: [Figure 1] 1 is a schematic diagram of an apparatus according to an embodiment. [Figure 2]1 is a flowchart of an element tracking method. [Figure 3] FIG. 1 is a schematic diagram illustrating an element tracking method according to an embodiment. [Figure 4] 1 is a flowchart of a density-assisted linking method according to an embodiment. [Figure 5] 1 is a flowchart illustrating a speed assistance linking method according to an embodiment. [Figure 6] 6(a) and 6(b) are diagrams illustrating the linking method of FIGS. 4 and 5. [Figure 7] 1 is a flowchart of a density-assisted linking method according to a further embodiment; [Figure 8] 8(a) to 8(d) are explanatory diagrams of the density-assisted linking method of FIG. 7. [Figure 9] 1 is a schematic diagram of an ultrasound imaging device according to an embodiment. [Figure 10] 10(a) to 10(e) are diagrams showing results obtained using the particle tracking method according to the embodiment. [Figure 11] 1 is a table of results obtained using particle tracking according to an embodiment. [Figure 12] 12(a) to 12(d) show results obtained using particle tracking methods according to embodiments. [Figure 13] 13(a) to 13(f) show results obtained using particle tracking methods according to an embodiment. [Figure 14] 1 is a table of results obtained using particle tracking according to an embodiment. DETAILED DESCRIPTION OF THE INVENTION

[0076] In general, super-resolution ultrasound imaging (SRUI) can be considered to consist of three processes: particle detection, particle localization, and linking. The following embodiments relate to a method for tracking elements, including particles, particularly microbubbles (MBs). Microbubbles are an example of a contrast agent, and MB detection can be achieved using many known methods, such as using particle probability images or comparing with a reference MB signal. Because MBs are much smaller than the imaging wavelength, they must be localized, such as by identifying the center of the MB signal. Then, the movement of the MB is tracked by linking the localized positions between consecutive frames.

[0077] FIG. 2 is a flowchart outlining the main steps of an element tracking method 10. This method 10 is intended to process medical images of an anatomical region of, for example, a human or animal subject, after administration of a suitable contrast agent, such as microbubbles, to the subject. The medical images are taken over a period of time, for example, a few seconds to a few minutes, so that the movement of the contrast agent through the anatomical region can be analyzed. Signals from the microbubbles can be detected using ultrasound techniques.

[0078] The imaging elements may be microbubbles or any suitable imaging agent, such as nanoparticles or contrast agent particles (CAPs). In other embodiments, suitable imaging techniques other than ultrasound may be used, such as CT scans, magnetic resonance (MR) scans, and positron emission tomography (PET) scans.

[0079] The first step 12 of the method 10 is acquiring a plurality of consecutive frames acquired using ultrasound imaging, each frame representing ultrasound image data.

[0080] The successive frames and image data represent an anatomical region of, for example, a human or animal subject, imaged over a period of time, such that each frame contains ultrasound image data representing the anatomical region at a different time.

[0081] Successive frames are captured by performing ultrasound measurements on a living human or animal subject. The ultrasound measurements indicate the presence of microbubbles administered to the subject. Data capture may occur simultaneously with frame processing, or data capture and frame processing may occur at different times. Successive frames may be stored and retrieved later for frame processing. An example of an ultrasound device for capturing successive frames is shown in Figure 9.

[0082] Microbubbles are contrast-enhancing agents that act as targets in ultrasound imaging. Microbubbles are typically gas bubbles less than 1 millimeter in diameter but larger than 1 micrometer. Microbubble solutions contain microbubbles of different sizes and shapes. Microbubbles typically range in diameter from 1 to 10 micrometers. Microbubbles are introduced into the subject using either a bolus injection or continuous infusion before and / or during the data acquisition phase.

[0083] In some embodiments, microbubbles are injected into a subject at a rate that allows for sparse distribution of particles within the frame. The appropriate injection rate may be experimentally determined or may be predetermined. The injection rate depends on many factors, including, for example, human physiology, microbubble suspension density, and imaging system. In some embodiments, the suspension density must be such that the microbubbles can be separated.

[0084] In some embodiments, the region being imaged is a portion of a subject's vascular bed. Microbubbles travel through the vascular bed, and imaging the microbubbles reveals the structure of the vascular bed. Typical vessels range in size from millimeters to a few microns.

[0085] In some embodiments, the vasculature comprises a blood vessel. The method uses tracking of the movement of microbubbles through the region to track the passage of blood to, from, or through an anatomical feature of interest. The anatomical feature may be a tumor or an organ.

[0086] Each frame is made up of pixels. Pixel size is typically 100 micrometers. In some embodiments, pixel size ranges from about 10 to about 1000 micrometers. The pixel size is therefore larger than the microbubbles being imaged. However, the signal from a microbubble typically spans the size of several pixels in the image. This is because it occupies the size of the point spread function (PSF), which is the response of the imaging system to the imaged microbubbles.

[0087] Each frame represents a view of a region of the subject. In a given view, for example, microbubbles may be seen overlapping in the resulting ultrasound image. As described elsewhere, the positional resolution can be selected to be greater than the pixel resolution.

[0088] Consecutive frames can be characterized by the duration they represent. Any duration suitable for collecting data to image the region of interest can be used. Additionally or alternatively, consecutive frames can be characterized by a frame rate. The frame rate can have any suitable value.

[0089] The second step 14 is directed to pre-processing the obtained successive frames. The second step 14 has two main components. The first component is an image registration process. This process is performed for each frame of the successive frames. In some embodiments, the image registration process is a rigid image registration. The image registration process operates to generate substantially aligned successive frames or video loops. The registration process operates to substantially remove image deformations or image artifacts due to externally induced motion, such as movement of the operator's probe. In other embodiments, alternative image registration processes are performed, such as non-rigid image registration. In some embodiments, the image registration process may be optional.

[0090] The second component of the preprocessing step is a filtering process performed on each successive frame. This process removes image artifacts, noise, and other speckles. The filtering process can be optional or partial, depending on the quality of the data and the image.

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

[0092] In some embodiments, the training process is optional, for example, the method may be performed using a preselected set of parameters without a training process.

[0093] Although shown generally as a separate method step, in some embodiments the training process forms part of other steps of the method.

[0094] The training process determines an optimized parameter set. The parameter set determined and used is described in more detail below. An initial manual evaluation is required to obtain the optimized parameter set. The training process may be performed on a subset of consecutive frames or on a separately provided training set of frames. The training data may be obtained separately.

[0095] Linking begins at step 18. Method 10 includes selecting a first frame from the sequence of frames and identifying a plurality of elemental signal portions within each frame at step 20. Step 20 corresponds to identifying one or more signal portions of the ultrasound imaging data of the frame as representing a microbubble or multiple microbubbles. Step 20 may also include classifying one or more identified signal portions as representing a single microbubble (single-microbubble signal portion) or multiple microbubbles (multiple-microbubble signal portions).

[0096] In step 22, location data is assigned to each component signal portion. Thus, location data is assigned to each single-microbubble signal portion and each multiple-microbubble signal portion. In this embodiment, assigning location data to the component signal portions includes fitting a mathematical function that determines location using at least one of the intensity, shape, and magnitude of the signal portions. In some embodiments, the mathematical function is a Gaussian function, and the determined location corresponds to a peak of the Gaussian function.

[0097] Information about the signal portions may be stored in a memory resource for later use by the link 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 classification, particle position (e.g., localization), particle position updates (e.g., via a Kalman filter), particle path (e.g., density map), particle velocity (e.g., velocity map), particle motion direction (e.g., one, two, or any suitable number of rose diagrams).

[0098] Following the identification and assignment of position data for the first frame, the process returns to step 18 to obtain a second frame. Element signal portion identification and position data assignment are then performed for the second frame. Steps 18 and 20 are repeated until all frames in the series, or a predetermined number of frames in the series, have been subjected to the detection and classification process.

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

[0100] In a subsequent step 24, signal portions representing microbubbles across different frames are linked using a linking method. The linking model takes as input at least a portion of the corresponding stored data for the frames. The linking model uses at least the 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. Linking signal portions between frames tracks the movement of microbubbles through anatomical regions of a human or animal subject. Frame linking is described in more detail below.

[0101] It will be appreciated that steps 18 through 22 are looped over all frames in the sequence. In some embodiments, steps 18 and at least part of step 24 are looped over a set of adjacent frames, e.g., a frame in the sequence and its immediate successor, such that the identification of signal portions (step 20), assignment of location data to each portion (step 22), and identification of links are looped (step 24, in part). More specifically, consecutive frames are ordered in time (1, ..., N), and steps are performed on a first pair of adjacent frames (1 and 2) to identify signal portions, assign location data, and obtain links between those frames, and then those steps are repeated on a second pair of frames (2 and 3) until all adjacent pairs have been processed.

[0102] FIG. 3 illustrates an element linking method according to an embodiment. It will be appreciated that in some embodiments, one or more steps, optionally all of the steps of FIG. 2, are included in the method of FIG. 3. For example, one or more steps, optionally all of the steps of FIG. 2, may be performed during each execution of FIG. 3. In other embodiments, the method of FIG. 3 is performed after an alternative particle detection and localization method. In FIG. 3, the particle detection followed by localization method is represented by step 32.

[0103] In the method of Figure 3, MB density is measured and criteria that each link must meet to be accepted are defined. These criteria form part of a density-assisted linking method, which can prevent links from being created between adjacent vessels. This is described in more detail with reference to Figure 4. Second, velocity and direction information within a local neighborhood is aggregated to guide the linking. This is called velocity-assisted linking and is useful for determining links between frames. In particular, the velocity-assisted linking procedure offers advantages in situations where a motion model is not available, such as when forming the first link of a track. This is described in more detail with reference to Figure 5.

[0104] The method of Figure 3 is implemented by running the particle tracking process twice. In step 33, an MB density map is obtained for use during the first and / or second run. The MB density map is described in more detail with reference to Figure 4.

[0105] In the first run (R1, 38), MB tracking is performed in stage 38, for example using a motion model. As will be explained in more detail with reference to FIG. 4, the linking process includes stage 38 using a motion model or nearest neighbor model to obtain potential links between frames. In stage 40, density-assisted linking is applied to refine the tracks and reject potential links that do not meet a density criterion. The result is a count of the tracks formed by the identified links.

[0106] In the second run (R2, 36), shown in the dashed box on the right, a velocity-assisted linking method is performed. In this embodiment, the linking process uses velocity information collected in Run 1. This may further refine the linking after density-assisted linking is applied, especially for the first link in a track. As described in more detail with reference to FIG. 5, the velocity-assisted linking method includes step 44, which uses a motion model (or other model) to identify candidate links between frames, followed by step 46, which uses a velocity-based or continuity criterion to refine the candidate links. In step 48, a maximum density search method is performed (step 52), which re-evaluates links rejected during the density-assisted linking of Run 1 to find whether alternative curved links exist. The maximum density search method is described in more detail with reference to FIG. 7. The result of Run 2 is final track information (step 52), including straight and curved tracks.

[0107] While FIG. 3 illustrates a linking method having a first run involving a first linking method followed by a second run involving a second linking method, where the results of the first run can be used as input for the second run, it will be understood that the linking methods described below (density seek method, velocity seek method, and maximum density seek method) can be performed independently or in a different order. In further embodiments, one or more linking methods can be repeated, and information from previous runs can be used in subsequent runs. For example, previously acquired knowledge about blood movement within blood vessels can be used. The process can be refined over multiple runs using different linking methods.

[0108] 4 illustrates a density-assisted linking method according to an embodiment. As will be explained below, this method links element signal portions in dependence on measurements of element density.

[0109] In steps 60 and 62, a number of link candidates are obtained using the location data of the microbubbles within the frame. The link candidates are formed in a known manner. Non-limiting examples of forming candidate links are given below. Candidate links may also be referred to as potential links.

[0110] In step 60, for each identified MB in a frame, a location in the next frame is predicted. In this way, for multiple identified MBs in the first frame, multiple predicted locations in the next frame are obtained. In this embodiment, the predicted locations are determined using a motion model or a motion continuity model. The identified MBs in frame t are labeled mi, where i=1, ..., M, where M is the number of identified MBs in frame t.

[0111] In step 62, the predicted locations of multiple MBs are compared with the locations of MBs detected in the next frame. In the next frame (t+1), the detected locations are labeled nj, where j=1, ..., N, where N is the number of MBs identified in frame t+1. A comparison is made between the predicted locations of MBs in the next frame and the identified locations of MBs in the next frame. In this embodiment, a cost matrix is ​​constructed and populated with a measure of the distance between each MB's predicted location and its identified location in the next frame. If M MBs are identified in the first frame and N MBs are identified in the second frame, the cost matrix is ​​an MxN matrix.

number

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

number

[0113] In the above question, mi'=(mx',my') and nj'=(nx',ny') are the coordinates of MBs. Link candidates or potential links are determined by minimizing each row of the cost matrix such that each MB is linked to and linked by one MB in two consecutive frames.

[0114] In step 64, a density map is obtained. The density map represents the density at multiple locations, e.g., multiple pixels. In this embodiment, the density map represents a measure of density across the entire frame for each pixel, in this embodiment, an average density value. To obtain the map, an average MB density value is calculated for each pixel and stored in a density map (also called a pixel density map). The average MB density is the number of MBs localized at this pixel across the entire image frame. Subsequent density criteria can use the density map for the entire frame or a subset, or alternative measures of density. In some embodiments, the density map may be obtained at another stage in the method or may be determined before step 60 of FIG. 4.

[0115] In step 66, density-assisted link criteria (also referred to simply as density criteria) are applied to the link candidates. In this embodiment, the criteria are applied to all possible links represented in the cost matrix. These criteria prevent links from being created between different vessels. To apply the density criteria, measurements of microbubble (MB) density are determined for multiple locations. In this embodiment, a measure of MB density is calculated for each of multiple pixels and for each frame of consecutive frames.

[0116] Density-assisted links can prevent links from being established between different vessels. For an image sequence under investigation, the MB density ρ(x,y) of a pixel (x,y) is defined as the number of MBs localized at this pixel across all image frames. As MBs move within vessels, correct links (mi and nj) between subsequent frames must satisfy a minimum MB density criterion, which updates the cost matrix. This cost matrix update rejects links that do not satisfy the density criterion.

[0117] The first density criterion uses the average element density across pixels along paths between pixels or other signal regions. In this embodiment, an average MB density is calculated along the path of the candidate link. Such a measure is calculated by averaging across MB density values ​​along the path. In this embodiment, these paths are straight paths. The first density criterion is applied by calculating the average MB density for each candidate link and comparing the average MB density to a minimum value. If the average MB density is below the minimum value, the link is rejected in step 68. In this embodiment, rejecting the link includes updating the cost matrix by assigning it a value of infinity.

[0118] In this embodiment, the first density criterion is expressed as follows:

number

[0119] The second density criterion used a measure of MB density in the form of the standard deviation of MB density across all pixels along the linear path. In this embodiment, the second density criterion is applied by calculating the average standard deviation for each candidate link and comparing the average standard deviation to a maximum value of the standard deviations. If the average standard deviation exceeds the maximum value, the link is rejected in step 68. In this embodiment, rejecting the link includes updating the cost matrix by assigning it a value of infinity.

[0120] In this embodiment, the second density criterion is expressed as follows:

number

[0121] The first criterion restricts links between MBs to meet a minimum MB density criterion. The second criterion ensures that MB density does not vary significantly along the trajectory. MB density measurements can be used to prevent or reduce links between different vessels, e.g., adjacent vessels. The density criterion ensures that there are a sufficient number of MBs and a smooth MB distribution along the trajectory, thereby avoiding MB jumps between vessels and improving vascular reconstruction.

[0122] In step 69, a number of tracks are determined using the identified links, each track being understood to be formed from a number of links. Once formed, the tracks represent the movement of the respective elements through a vessel, for example, following blood flow.

[0123] Figure 5 illustrates a velocity-assisted linking method according to an embodiment. As will be explained below, this method links constituent signal portions in dependence on a velocity measurement.

[0124] In steps 70 and 72, a motion model is used to identify possible links between frames. Steps 70 and 72 correspond to steps 60 and 62 in FIG.

[0125] In step 74, a velocity map and a direction map are obtained. In this embodiment, the velocity and direction map are obtained using previously determined track information. For example, as depicted in FIG. 3, the velocity and direction map is obtained during a second run. This second run occurs following a first run in which trajectory information was obtained, and this trajectory information is used to determine the velocity and direction map. All determined MB tracks are used to generate the velocity and direction map. In other embodiments, the tracks may use alternative methods other than density-assisted linking. In general, although FIG. 3 depicts a first and second run, either or both runs may be performed iteratively using information obtained from previous runs. In this manner, the resulting tracks may be refined over multiple iterations.

[0126] The velocity map represents the average velocity of the previously determined tracks at each of a plurality of locations, e.g., each of a plurality of pixels. In this embodiment, the average velocity of all tracks is calculated for each pixel. The average velocity pixel values ​​form the velocity map. In addition to the average velocity map, a map of the standard deviation of the velocities of all tracks is also generated at this stage. The velocity standard deviation map is determined by calculating the standard deviation σV(x,y) of the velocities of all tracks passing each pixel.

[0127] In step 74, a direction map is also determined using the track information. The direction map represents the average direction of the previously determined tracks at each of a plurality of locations, e.g., at each of a plurality of pixels. In this embodiment, the average direction of all tracks is calculated for each pixel. The pixel value of the average direction forms the direction map. In addition to the map of average directions, a map of the standard deviation of the directions of all tracks is also generated at this stage. The map of standard deviation of the directions is determined by calculating the standard deviation of the directions of all tracks that pass through each pixel σθ(x,y).

[0128] In step 76, a velocity-based criterion is applied to the candidate links identified in step 72. The velocity-based criterion is also called a continuity criterion. As MBs move along the vasculature, their velocity and direction change between adjacent frames, but should be within a certain range of continuity. In this embodiment, two velocity-subcriteria are applied: fixed continuity and variable continuity. The criterion that yields the best results may depend on the results obtained in R1. The fixed continuity criterion may provide better results in situations where all tracks across the image have fairly good continuity in both velocity and direction, with small and well-defined local and global variations. In some embodiments, the application of the continuity criterion may represent the rejection of physically improbable or impossible microbubble paths.

[0129] In this embodiment, first criteria related to link speed and direction are applied. Applying the speed criteria involves first comparing the link speed with a lower limit (corresponding to the minimum speed) and an upper limit (corresponding to the maximum speed). Applying the direction criteria involves comparing the link direction with a maximum angle relative to a predetermined direction and a minimum angle relative to a predetermined direction. If the link speed is not within the range defined by the maximum and minimum speeds, or if the direction is not within the range defined by the maximum and minimum angles, the link is rejected.

[0130] The first criterion can be expressed as follows:

number

[0131] In this embodiment, a second criterion is applied regarding the link speed and direction. The second criterion is applied when the continuity of speed and direction is inconsistent. In this embodiment, the variation criterion is applied based on a previously determined measure of continuity in the track. In this embodiment, the standard deviation of speed and direction is determined for each track, for each pixel, and used to apply the criterion. The second criterion includes comparing the measured link speed with upper and lower bounds determined from the per-pixel speed and the corresponding per-pixel standard deviation.

[0132] The second criterion is expressed as follows:

number

[0133] As described above, a link is rejected because it does not satisfy any of the criteria. In this embodiment, rejecting the link includes updating the cost matrix by assigning it a value of infinity. In step 80, tracks are formed using the cost matrix, e.g., as described with reference to step 69 of Figure 4.

[0134] Figure 6(a) illustrates the concept of density-assisted links. Figure 6(b) illustrates the concept of velocity-assisted links. Figures 6(a) and 6(b) together show a vasculature 102 having a first vessel 104 and a second vessel 106. The gray arrows indicate the direction of blood flow through the first and second vessels. Figures 6(a) and 6(b) show microbubbles in different frames.

[0135] FIG. 6(a) shows MB 108 in a first frame (time t) and four MBs in a second subsequent frame (time t+1). The four MBs in 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 between each MB in the second frame and the MB in the first frame are depicted in FIG. 6(a). These may be determined, for example, using a motion-based model or nearest neighbor model, as described above. Four such potential links (112a, 112b, 112c, and 112d) are depicted in FIG. 6(a).

[0136] Similarly, Figure 6(b) shows the first MB 114 in the first frame (time t) and four MBs in the second frame (time t+1). The four MBs in the second frame are labeled 116a, 116b, 116c, and 116d. The velocity-assisted linking method limits the number of potential links to a single potential link 118 between MB 114 in the first frame and third MB 116c in 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 ring. The sector has a subtraction angle defined by the allowable angle range described above. The thickness of the sector (between the inner and outer radii) is defined by the lower and upper velocity limits. It can be seen that this shape allows a physically likely link (e.g., link 118) to be selected from the candidate links.

[0137] FIG. 7 is a flowchart of a maximum density method according to an embodiment. As described with reference to FIG. 4, a linking method, such as a density-assisted linking method, can reject a straight link between locations. If a straight link is accepted, it can be used to form part of a track. If a straight link is rejected, it may be incorrect. Alternatively, the link may have the correct start and end points, but may not be straight, but may have a non-linear trajectory. The maximum density-based linking method is described below. It will be understood that the method of FIG. 7 is performed with reference to a pixel map, as described above.

[0138] In step 82, the start and end positions of the rejected links are obtained. In this embodiment, these correspond to the initial pixel and the final pixel, labeled mi and nj. In step 84, a neighborhood window is applied to the initial pixel. The neighborhood is a 3x3 array of pixels, which defines eight neighborhoods.

[0139] In step 86, one of the eight neighboring pixels is selected based on at least the MB density value, e.g., obtained from a pixel map. In this embodiment, the neighboring pixels are filtered based on a threshold density value, and only pixels equal to or greater than the threshold are considered. The filtered pixel with the closest distance to the end of the link is then selected as the next pixel on the new path.

[0140] In step 88, a new neighborhood is defined around the pixel selected in step 86, and further pixels are selected. At this stage, the window is slid to a new position, and eight new neighboring pixels are defined around the selected pixel. Selection of one of the neighboring pixels is performed similarly to step 88, by filtering the pixels and selecting the filtered pixel with the closest distance to the end point. The further selected pixel is added as the next pixel on the new path.

[0141] In step 90, the selection of pixels on the new path is repeated until a final pixel is reached. When the final pixel is reached, a new path is defined through the set of selected pixels. The new path has a non-linear trajectory. Each pixel is labeled with the sequence (x1, y1)...(xN, yN). If (xN, yN) does not reach the terminal pixel (nj) in this way, the method ends and the link is rejected as there is no alternative path.

[0142] Figure 8 illustrates the steps of the maximum density method. Figure 8(a) shows a pixel map with shading representing the density of each pixel value. For Figures 8(a) through 8(d), pixel density is shown on a scale of nine possible values, ranging from zero density (first density 202) to the 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, eighth density 216, and ninth density 218. While Figure 8 depicts the densities on a scale, it will be understood that these labels are provided for illustrative purposes only and that densities may be implemented on a continuous scale. The same density labels as in Figure 8(a) also apply to Figures 8(b) through 8(d) but have been omitted for clarity. The pixels can be divided into high density pixels (corresponding to labels 210, 212, 214, 216, 218), low density pixels (corresponding to labels 204, 206, 208), and zero density pixels (corresponding to label 202). Zero density means that virtually no MBs are detected.

[0143] The filled circle 220 represents the start pixel and the filled triangle 222 represents the end pixel. The filled circle corresponds to mi and the filled triangle corresponds to nj. Figure 8(a) shows the first link candidate, a straight link 224. As part of the maximum density method, this candidate link is rejected because it violates the density criterion. This straight link lies between mi and nj and is rejected by the path density criterion.

[0144] 8(b) illustrates the examination of the neighborhood of the starting pixel 220 (corresponding to step 84). The neighborhood is a 3x3 pixel region 230 centered on the starting pixel 222. Of the eight neighboring pixels in pixel region 230, the pixel closest to the ending pixel 220 (greater than or equal to a predetermined threshold) is identified as the first identified pixel 228 of the new path and assigned coordinates (x1, y1).

[0145] 8(c), the process is repeated by examining a further neighborhood around the initially identified pixel 228, corresponding to step 86. The further neighborhood defines a further 3x3 pixel region 232 centered on the initially identified pixel 228. The pixel closest to the final pixel 220 that is greater than or equal to the pixel threshold is identified as the second identified pixel 234 and assigned coordinates (x2, y2).

[0146] This process is repeated until the following pixels are identified, as depicted in FIG. 8(d): 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, which are labeled (x1, y1)...(x6, y6). In this example, (x6, y6) equals final pixel 222, so the method ends. A curved link 244 is then fitted between start pixel 220 and final pixel 222 through the identified pixels to form a new link.

[0147] Determining the curved links between the predetermined initial and final portions can be viewed as finding a path based on density values ​​and proximity to rejected ends. Such a path may include adjusting a speed scale according to the determined path. For example, if a path is determined to be curved and therefore longer than a straight path, increasing the path's speed. A straight link may represent a physically impossible path in the sense that there are restrictions or other obstacles along the path. A curved link instead of a straight link may represent a physically possible path in the sense that there are no restrictions or obstacles along the path, or that the path is a path. For example, a rejected straight link may represent a path between different containers or a path that crosses a container boundary, while a curved link may represent a path contained within a single container.

[0148] One embodiment of an ultrasound imaging device 100 is shown schematically in FIG. 9. A linear array transducer 110 is configured to transmit energy to an imaging target (e.g., a portion of a human or animal body) and receive ultrasound echoes from the target. In other embodiments, any suitable transducer can be used to receive 2D or 3D ultrasound data and need not be a linear array. For example, the transducer may be a phased array, a curvilinear array, or other suitable array.

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

[0150] The digitized ultrasound data is processed by a processor 114, and the resulting processed data may also be stored in memory 116, or in a separate memory. The processor 114 may be configured to perform beamforming of the digitized echo data. In some embodiments, multiple processors 114 may be used.

[0151] The ultrasound imaging device 100 also includes a screen 118 for displaying the ultrasound image (which is the result of further post-processing to suit display requirements) and one or more user input devices 120 (e.g., a keyboard, mouse, or trackball) for receiving input from a user of the ultrasound imaging device 112, e.g., a sonographer.

[0152] 9, the ultrasound imaging device 100 is configured to acquire and process ultrasound data using a linear array transducer 110. In other embodiments, a separate processing device (e.g., a workstation or general-purpose computer) can be used to process ultrasound data previously acquired by the ultrasound imaging device.

[0153] In some embodiments, the organs imaged by the imaging system can include any human or animal organ, such as an organ of the digestive system, including the liver, pancreas, etc., an organ of the urinary system, including the kidney or bladder, an organ of the cardiovascular system, including the heart, any sensory organ, or an organ of the central nervous system, including the brain. Anatomical features can include any component of the lymphatic system, including lymphatic vessels and lymph nodes, for example. Anatomical features can include any component of the cardiovascular system, including the heart, arteries, veins, capillaries, or portions of a vascular bed, for example.

[0154] As microbubbles move within vasculature, a first application is to track the microbubbles to visualize the vasculature, as described above. However, vasculature also moves due to subject movement, e.g., heartbeat and breathing. Thus, the track contains information not only about blood movement, but also about the movement of the vasculature walls due to vasculature pulsatile motion propagated by the heart. Therefore, a further application is to use the trajectory to perform measurements or to represent the pulsatile motion of the subject.

[0155] In some embodiments, for example in CT, MR or PET applications, contrast agents, in contrast to microbubbles, may travel outside the vascular lumen, and therefore another application is to track the movement of contrast agents outside the vascular lumen, which may be slow.

[0156] The above-described embodiment describes the step of identifying signal portions. Further comments on non-limiting examples of identifying microbubbles or element signal portions are provided below. As a first example, a binary image is generated from the original grayscale image. In this step, multiscale Haar-like features that measure local contrast at different shapes and sizes are generated. The Haar-like features are formed using non-local and statistical mapping of the original grayscale image, leading to a binary image. Using the generated binary image, one or more signal portions that represent a microbubble or multiple microbubbles are identified in the generated binary image. Signal portions that represent microbubbles of all sizes and intensities are identified. The signal portions identified in this stage are candidate signal portions. Some candidate signal portions may correspond to noise. There is no limit to the number of signal portions that can be identified. Therefore, this stage identifies all candidate signal regions, regardless of their characteristics. The final number of signal regions, and therefore the number of represented microbubbles, is determined by subsequent filtering and selection stages. A first number of signal portions / microbubbles identified in a first frame can be different from a second number of signal portions / microbubbles identified in a second frame, and the first and second numbers can be independent.

[0157] The generated binary image is then subjected to a first filtering and / or thresholding process, which acts to filter out signal portions that do not correspond to microbubbles or meet other criteria. The first filtering process filters out isolated detected pixels. A strict or less strict connectivity filter is applied. The first filtering process also includes thresholding. The size of the signal region is determined and compared to a predetermined threshold. If the size of the signal region is smaller than a minimum size, the signal region is discarded. The intensity of the signal region is also determined and compared to a predetermined threshold. The intensity of the signal region can be determined, for example, by summing the intensities of all pixels in the signal region. If the intensity of the signal region is lower than a threshold intensity value, the signal region is discarded. Depending on the quality of the data and image, the filtering and / or thresholding may not be performed, or may only be performed partially.

[0158] Next, a particle probability image (PPI) is generated from the binary image. Each detected signal portion is enhanced. Each signal portion in the binary image is refined based on the foreground and background values ​​of the PPI. An image smoothing process is then applied to the original grayscale image. The smoothing process involves convolving the original image with a Gaussian smoothing kernel. Furthermore, local maxima are found and refined using the following criteria: the local maxima must be within the segmented region and have a PPI value greater than or equal to a predetermined particle region threshold. A second filtering process is applied. This involves generating a watershed transform and discarding detected signal portions that are too large. A location is then assigned to each remaining signal portion. A geometrically weighted centroid is determined for each signal portion using both the signal portion's size and intensity. The size, shape, and intensity values ​​may be determined in this step or reused from previous steps, such as step 36. The locations may be assigned at a higher resolution than the resolution of the original image. An optional step is to classify the identified signal portions as corresponding to single or multiple microbubble events. This classification is based on the size, intensity, and shape of the signal portion. Signal portions corresponding to multiple microbubbles are assigned only a single location, but may remain classified as multiple microbubble events to allow for multiple paths and more accurate particle density assignment in the later linking stage.

[0159] Furthermore, the above embodiments are described with respect to a motion link model. In the following description, a motion model is provided, but it will be understood that alternative links using other motion models or nearest neighbor models may be used. In alternative embodiments, a different motion model is implemented as part of the energy matrix update. 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.

[0160] The first stage of the linking step using the motion model involves linking identified single or multiple microbubble signal portions to track the movement of microbubbles between successive frames and forming track segments using position data of the single and multiple microbubble signal portions. The second stage links the formed track segments to generate multiple microbubble tracks. The first stage of linking signal portions between frames to form track segments is based on a parameter, maximum displacement (labeled MD), which controls the number of pixels a particle is allowed to move from one 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 neighbor method between signal portions of successive frames. Each neighborhood region is defined by a value, MD, which is the maximum allowable displacement between frames. The energy matrix can include information such as particle intensity, size, and shape to facilitate particle recognition at the next position. The linking model is operable to link at least one single microbubble signal portion in at least one of the frames to multiple microbubble signal portions in at least one other subsequent frame, or vice versa. The linking model is also operable to link multiple microbubble signal portions representing different numbers of microbubbles in different frames. After an initial linking stage, the energy matrix is ​​optimized and updated using multiple Kalman filters and the assumption that the microbubbles represented by the signal portions move according to a linear motion model.

[0161] We offer the following non-limiting comments on our experimental results: To test the new linking strategy described above (also called PTNS—Particle Tracking with Neighborhood Similarities), we applied the method to three different datasets: synthetic images, animal images, and human prostate CEUS images. These datasets allow us to quantitatively evaluate the performance of PTNS under a variety of different conditions, some of which have well-defined, known features that can be compared to the structural and mechanical features recovered by PTNS.

[0162] First, a synthetic flow model is created to simulate blood flow within a vascular network. Next, point-like particles are injected into the network and move along with the flow. These particles are blurred using a variable point spread function to mimic the MBs observed in real in vivo CEUS data. Finally, white Gaussian noise is added to each frame. This results in a realistic CEUS dataset of MBs moving at various speeds along the predesigned network structure. Each MB is assigned an ID number and used as ground truth for comparison with PTNS tracking results.

[0163] To evaluate the performance of PTNS, we compare the speed maps of MBs obtained by different particle tracking methods with the ground truth shown in Figure 10. The velocity map records the average velocity of all tracks passing through each pixel in the image sequence and thus displays not only the flow dynamics but also the structure of the vascular 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).

[0164] We apply the nearest neighbor linking method combined with a motion model to synthetic data. As shown in Figure 10(b), the overall structure of the network is nearly present, but some structural details are missing. The velocity in the center of the network is close to the ground truth, but significantly lower at both ends. This is because the nearest neighbor method tends to link nearby MBs. Even with a motion model, this method can result in large link errors, especially in high-speed environments.

[0165] PTNS is applied to synthetic data. The results of R1 are shown in Figure 10(c), which show a significant improvement in recovering both the network structure and flow dynamics compared to the ground truth. This must be attributed to the density-aided link, which is the main difference from the nearest neighbor method. However, erroneous links are still noticeable in R1, mainly in the central region of Figure 10(c), where the tracks are more densely packed. Since the results of R1 show that the speed continuity is nearly uniform throughout the network, we can apply the speed-aided link model (Equation 5) for a second run.

[0166] As shown in Figure 10(d), most of the incorrect links in the middle region of R1 have been removed because linking was performed based on the velocity information collected by R1. The maximum density search method was applied to include the bent links, and the results are shown in Figure 10(e). Visually, it does not appear to be significantly different from 10(d), but the quantitative analysis in Table 1 shows that the bent links improve the results.

[0167] Figure 11 shows a table (Table 1) showing the number of correct links (CL) generated by the linking method, the total number of links (DL) generated by the method, and the total number of links in the ground truth (GT). From these, two statistics, accuracy and Jaccard index, are calculated to quantitatively measure the performance of different linking methods. From left to right in the table, the methods are: nearest neighbor + motion model method, PTNS R1, PTNS R2 with straight links, PTNS R2 with bent links, and PTNS R2 with a combination of straight and bent links. It can be observed that applying PTNS improves both accuracy and Jaccard index. Applying the maximum density search method to R2 slightly decreases accuracy, but significantly increases the Jaccard index. Therefore, PTNS not only improves tracking accuracy but also allows for more efficient use of data, which is important in clinical settings.

[0168] PTNS was also tested on an animal dataset. Data were acquired from sheep ovaries using a bolus injection of 1.2 mL of contrast agent, SonoVue (Bracco, Geneva, Switzerland). 936 CEUS frames were collected at a frame rate of 5 Hz. Because this is an in vivo dataset, no ground truth exists. However, the simplicity of the data allows for a sufficient assessment of the performance of PTNS. It also contains very salient features that challenge the ability of PTNS to reconstruct structures using MB tracks. Therefore, PTNS was applied to this dataset as described with reference to Figure 3. The results are shown in Figure 12.

[0169] The MBs in each CEUS frame are located and collected into an MB number map, where each pixel (x, y) in the MB number contains the total number of MBs localized at that pixel for all image frames, as described previously in Section 2.2. The results are shown in Figure 12(a). Three distinct features are clearly visible in the MB number map. These are labeled 1, 2, and 3 in Figure 12(a) and are a horseshoe-shaped structure, a curved structure, and a thin vessel. Reproductions of these three features are used in place of the ground truth for this dataset. As with the synthetic data, speed maps can be used to display dynamics and structure, allowing for a clear visualization of the performance of each tracking method.

[0170] The nearest neighbor linking method is applied to the sheep data. Figure 12(b) clearly shows that feature 1 is not well reproduced. The horseshoe shape is heavily smudged, indicating that many links are made across the gap in the middle of the structure. Feature 2 is mostly well reproduced, but there is an error in the link connecting it to feature 1 on the left. Of the three features, feature 3 is reproduced best. This region has low MB concentration, slow velocity, and thin structure, conditions that nearest neighbor methods excel at.

[0171] Next, we apply PTNS to the sheep data. Figure 12(c) shows the speed map of R2 using only straight links, and Figure 12(d) shows the speed map of R2 using a combination of straight and bent links. Features 1 and 2 are better reproduced using PTNS than by nearest neighbor. In Figure 12(c), the blank area in the middle of the horseshoe shape of feature 1 is preserved, and the erroneous link connecting feature 1 to the left of feature 2 is eliminated. However, straight links alone cannot preserve the connection between feature 1 and feature 2 on the right side of feature 2. Combining straight and bent links reproduces this connection. Furthermore, the structure of features 1 and 2 is further improved, as seen in Figure 12(d).

[0172] In Figure 12(c), feature 3, the thin blood vessels, are not well reproduced. Figure 12(d) shows some improvement, but there are still missing parts. Because the MB density is sparse in this region, the density-assisted linking is likely to reject links. Because of this sparseness, there is no clear path with the highest MB density for the maximum density search method to follow, and as a result, it is unable to find alternative paths that are not linear. As a result, we can see that the maximum density search method is not suitable for datasets with very low MB concentration.

[0173] We next applied PTNS to human prostate data. CEUS data were collected with full ethical approval from patients scheduled for radical prostatectomy at the Western General Hospital in Edinburgh. Patients were injected with contrast at a near-constant flow rate, and data were collected for 3–4 minutes at a frame rate of 10 Hz using an iU22 Philips scanner with a C10-3V transducer. As with the sheep dataset, no ground truth exists, but specific structural and mechanical features and post-pathological examination results can be used to demonstrate the performance of the tracking method. The results of this test are shown in Figure 13.

[0174] Figure 13(a) shows a B-mode prostate image of the patient under investigation, with the prostate boundary outlined in yellow. Figure 13(b) shows a MB count map created by collecting all MBs from 2,070 CEUS frames of the prostate. The MB count density pattern in the map shows a vascular network, though without any dynamic information. PTNS tracks these MBs sequentially, starting with R1. Figure 13(c) shows the track count map (the number of tracks passing each pixel) of R1, and Figure 13(d) shows the corresponding velocity map. The cancerous region was identified by the pathology team and is enclosed by a red dashed line in Figure 13(c). As can be seen, the cancerous region exhibits a high track count and high velocity. Even in healthy regions, the high track count and high velocity are evident, particularly in the central region. However, the cancerous region does not exhibit a distinct spatial structure compared to healthy regions.

[0175] Next, PTNS R2 is applied with a variable range of continuity, as defined in Equation 6. The results of different continuity levels were investigated by adjusting the control parameter y. Figures 13(e) and 13(f) show the track count and velocity maps for R2 with y = 1.5. Comparing Figures 13(e) and 13(c) reveals that R2 reduces the total number of tracks in the prostate, and the structure resulting from the imposition of the continuity constraint is also clearly visible. Similar features due to the continuity constraint can also be observed by comparing Figures 13(f) and 13(d).

[0176] It is known from the literature that the vasculature of prostate cancer differs from that of healthy tissue (Forster et al., 2017; Alizadeh et al., 2013; Vaupel and Kelleher, 2012). Cancer vessels are generally larger than healthy vessels (Forster et al., 2017; Muller et al., 2008), vessels are increasingly tortuous (Alizadeh et al., 2013), and blood flow is heterogeneous (Vaupel and Kelleher, 2012; Jochumsen et al., 2020). Therefore, we expect that varying the continuity range will result in slight differences in structure and dynamics between cancer and healthy regions. To quantify this, we compared the number of tracks in cancer and healthy regions at different levels of continuity restriction.

[0177] Figure 14 contains a table (Table 2) comparing the number of tracks in the cancer and healthy regions of a patient (SRI010) with respect to two different y-views, R1 and R2. As can be seen in Table 2, both the cancer and healthy regions show a decrease in the number of tracks as the range of continuity becomes more restrictive, but the decrease in the number of tracks is larger in the cancer region: 5% at y=1.5 and 9% at y=1.

[0178] At y=1.5, the decrease in track count is 5% and at y=1, it is 9%. Thus, the cancerous region exhibits a larger decrease in track count than the healthy region, implying that the former exhibits less continuous motion than the latter. This indicates that the blood vessels and corresponding blood flow in the cancerous region behave more chaotically in structure and dynamics, respectively. This finding is consistent with previous observations on prostate cancer vasculature (Alizadeh et al., 2013; Vaupel and Kelleher, 2012; Jochumsen et al.). These findings indicate that PTNS may be a useful tool to identify distinctive structural and mechanical features of cancerous regions to aid in the diagnosis of prostate cancer.

[0179] A novel method, particle tracking using neighborhood similarity, attempts to address several key issues with current particle linking approaches and methods when applied to CEUS images. Using synthetic data, PTNS has been shown to improve particle tracking performance. By improving accuracy and Jacqard exponent, PTNS uses data more efficiently than previous methods. Using animal data, PTNS demonstrated its ability to reconstruct complex structures in in vivo data, despite the lack of ground truth. Finally, PTNS was applied to human prostate data to investigate the structure and dynamics of the vasculature in patients with prostate cancer. Because PTNS is a two-run process, different levels of motion continuity can be explored in the second run by varying the control parameters. Using this method, we observed that cancerous regions exhibited more chaotic vascular structures, and consequently, a more restrictive velocity continuum resulted in a greater reduction in track counts. This is consistent with our current understanding of prostate cancer. These results suggest that super-resolution ultrasound imaging of prostate cancer has the potential to be developed into a diagnostic and localized treatment tool.

[0180] Although specific embodiments have been described above, it should be understood that these embodiments are merely illustrative and that the claims are not limited to these embodiments. Those skilled in the art will be able to make modifications and alternatives to the described embodiments that are contemplated as falling within the scope of the appended claims. Each feature disclosed or illustrated herein can be incorporated into any embodiment, either alone or in any suitable combination with other features disclosed or illustrated herein. In particular, those skilled 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 provide advantages or benefits when used separately from one or more of the other features of the embodiments of the present disclosure, and that different combinations of features are possible other than the specific combinations of features of the embodiments of the present disclosure described above.

Claims

1. 1. A method for tracking an element, comprising: obtaining, for each successive frame, position data including a respective position assigned to each of a plurality of elementary signal portions within a frame; using the acquired position data to link element signal portions in one frame to element signal portions in other frames to track element movement through an area of ​​a subject; Including, a) the linking of said elementary signal portions depends on measurements of elementary density in at least a portion of said one frame and in at least a portion of said other frame; or b) the linking of said element signal portions depends on measurements of element velocities in at least a portion of said one frame and in at least a portion of said other frame; A method characterized by:

2. 2. The method of claim 1, wherein said linking of elementary signal portions forms a plurality of links, each of said plurality of links linking a respective pair of elementary signal portions.

3. obtaining a plurality of potential or candidate links between a first frame and a second frame of the consecutive frames; selecting one or more of the potential or candidate links depending on element density measurements and / or element velocity measurements of the first and / or second frames; The method of claim 1 or 2, further comprising:

4. obtaining the plurality of potential or candidate links comprises applying a motion model and / or a nearest neighbor model and / or a known link model; and / or 4. The method of claim 1, wherein the step of selecting one or more of the potential or candidate links comprises applying criteria to the potential or candidate links to obtain physically possible and / or likely links between frames.

5. the links of the elementary signal portions form a plurality of tracks; Each of the plurality of tracks may include a respective plurality of links; 5. A method according to any preceding claim, optionally wherein each of said plurality of tracks is a representation of the movement of a respective element.

6. the elements include particles and / or contrast elements, e.g., microbubbles; 6. The method of any one of claims 1 to 5, wherein the microbubbles optionally comprise bubbles having a diameter of at least 1 micrometer but less than 1 millimeter, and optionally at least 1 micrometer but less than 10 micrometers.

7. acquiring successive frames containing ultrasound or other medical image data representing anatomical regions of a human or animal subject at respective different times; for each of said frames, identifying a plurality of elementary signal portions and assigning respective position data to each of said elementary signal portions; The method of any one of claims 1 to 6, further comprising:

8. the linking of the element signal portions is dependent on a measure of element density; 8. The method of claim 1, further comprising the step of optionally determining a measure of element density for each of a plurality of locations, e.g., for each of a plurality of pixels, and / or for each of the successive frames.

9. The step of determining a measure of element density comprises: a) using a density map representing density at multiple locations, e.g., multiple pixels; or b) using the density for all frames of said consecutive frames; or c) using density for a subset of said consecutive frames; 9. The method according to claim 7 or 8, comprising:

10. The element density measurement is density per pixel, average element density, the average element density over pixels along the path between element signal regions, the standard deviation of the element density and / or the standard deviation of the element density over pixels along the path between element signal regions 10. The method according to claim 1, further comprising:

11. the linking step includes a density-assisted linking method; The density-assisted link method includes: applying at least one density criterion to the potential links, e.g., a minimum density criterion; 11. The method of claim 1, further comprising the step of: optionally rejecting potential links that do not meet at least one density criterion.

12. 12. The method of claim 11, wherein the at least one density criterion includes a minimum value of average element density along a path, e.g., a straight path, and optionally rejecting a link if the average element density along the path of that link is below the minimum value.

13. 13. The method of claim 11 or 12, wherein the at least one density criterion comprises a maximum value of the standard deviation of the element density along a path, e.g. a straight path, and optionally a link is rejected if the standard deviation of the element density along the path for that link exceeds the maximum value.

14. 14. A method according to any preceding claim, wherein the element velocity measurements comprise element velocity measurements and / or element orientation measurements.

15. determining a measure of element velocity; 15. A method according to any preceding claim, wherein optionally a measure of element velocity is determined for each of a plurality of positions, e.g. for each of a plurality of pixels, and / or for each frame of successive frames.

16. 16. The method of claim 1, wherein the linking step comprises a speed-assisted linking method.

17. 17. The method of claim 16, wherein the speed-assisted link method includes obtaining a speed map and / or a direction map.

18. a) the speed map and / or direction map are generated using a previously determined track, or b) the map of velocities represents an average velocity of a previously determined track at each of a plurality of locations, e.g., at each of a plurality of pixels; or c) the map of orientations represents the average orientation of the previously determined tracks at each of a plurality of locations, e.g., at each of a plurality of pixels; 18. The method of claim 17.

19. The speed assistance link method includes: applying at least one continuity criterion to the potential links; 19. A method according to any of claims 16 to 18, optionally rejecting potential links that do not satisfy at least one continuity criterion.

20. 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 to a predetermined direction, or a minimum angle to a predetermined direction, and / or depends on a previously determined measure of continuity in the track, such as a standard deviation of the speed and / or a standard deviation of the direction.

21. applying the at least one continuity criterion to the potential links, defining a partial sector or other shape from the elements based at least on velocity and / or angle; determining whether the potential link resides at least partially in a partial sector or other shape; 21. The method of claim 19 or 20, comprising:

22. 22. A method according to any preceding claim, wherein the linking step comprises using velocity and / or direction for a subset of consecutive frames and / or for all frames.

23. the linking step includes determining a maximum density; Optionally, determining the maximum density comprises determining a curvilinear link between two elementary signal portions; 23. A method according to any one of claims 1 to 22, further optionally wherein the step of determining the curvilinear link comprises the step of finding a path with a maximum element density between two element signal portions in two consecutive frames.

24. determining the maximum density includes using a map of densities to find a path of maximum element density; 24. A method according to claim 23, wherein optionally the density map comprises a respective element density for each of a plurality of locations, e.g. for each of a plurality of pixels.

25. determining the maximum density includes adjusting a scale of velocity along the determined path; 25. A method according to claim 23 or 24, characterized in that if the path is determined to be a curved path and longer than a straight path, the speed along the path is increased.

26. 26. A method according to any preceding claim, wherein the linking step comprises replacing physically impossible or improbable links, e.g. straight links, with physically possible or possible links, e.g. non-straight links.

27. The linking step includes: determining a set of links using nearest neighbor and / or motion model methods; using element density measurements and / or element velocity measurements to refine the set of links; 27. The method of claim 1, comprising:

28. at least some of the elements are present in the vasculature of a human or animal; 28. The method of any preceding claim, further comprising the step of tracking the path of at least some of the vessels by tracking the movement of the elements through the region.

29. The links of the elementary signal parts form a plurality of tracks, 29. A method according to any preceding claim, wherein each of the tracks comprises a respective plurality of links.

30. 1. An image processing system, comprising: obtaining, for each successive frame, position data including a respective position assigned to each of a plurality of element signal portions within a frame; using the acquired position data to link element signal portions in one frame to element signal portions in other frames to track movement of elements through an area of ​​a subject; a processing resource configured to: a) the linking of said elementary signal portions depends on measurements of elementary density in at least a portion of said one frame and in at least a portion of said other frame; or b) an image processing system, wherein the linking of the element signal portions depends on measurements of element velocities in at least a portion of the one frame and in at least a portion of the other frame.

31. 1. An imaging system comprising: an ultrasound or other scanner configured to perform a scan of a human or animal subject to acquire successive frames; 31. An image processing system according to claim 30 configured to receive and process the successive frames to track movement of an imaging element through a region of the subject; 1. An imaging system comprising:

32. 30. A computer program product comprising computer readable instructions capable of carrying out the method of any of claims 1 to 29.

Images