Accelerated coupled filtering method and system for tissue deformation analysis

Abstract

There is provided an accelerated coupled filtering method for tissue deformation analysis. The method includes initializing, by a first processing device, a plurality of first motion parameters; applying, by a second processing device, a first filter on a pre-deformation image of a tissue to obtain a filtered pre-deformation image, wherein the second processing device is different from the first processing device; applying, by the second processing device, a second filter on a post-deformation image of the tissue to obtain a filtered post-deformation image, wherein the filtered pre-deformation image and the filtered post-deformation image are correlated by a first motion matrix comprising the plurality of first motion parameters, and wherein the plurality of first motion parameters comprise at least three fundamental first motion parameters, each of the three fundamental first motion parameters represents movement of the tissue along an axial direction relative to the axial direction, an elevational direction and a lateral direction, respectively; estimating, by the second processing device, a respective value for each of the plurality of first motion parameters, wherein each of the estimated respective values represents a difference between the pre-deformation and post-deformation filtered images; and in response to determining that at least one of the estimated respective values meets at least one predefined criterion, updating the at least one of the estimated respective values as an optimal value for at least one corresponding first motion parameter of the plurality of first motion parameters.

Description

ACCELERATED COUPLED FILTERING METHOD AND SYSTEM FOR TISSUE DEFORMATION ANALYSISCROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application claims priority from U.S. Provisional Patent Application No. 63 / 666,678, filed on July 2, 2024, which is hereby incorporated by reference in its entirety.FIELD OF INVENTION

[0002] The present invention relates broadly, but not exclusively, to accelerated coupled filtering methods and systems for tissue deformation analysis.BACKGROUND

[0003] Ultrasound image-based tissue deformation analysis can be used to measure tissue stiffness and reveal useful information for clinical diagnosis. Such analysis is typically performed by comparing images taken before and after deformation of tissues. However, feature motion decorrelation caused by the ultrasound imaging procedure can greatly impede accuracy. To address this problem, a coupled filtering method was proposed and implemented to compensate for feature motion decorrelation analytically. Although the coupled filtering method can achieve much higher accuracy as compared to other existing methods, such method is computationally expensive. Thus, the implementation of such method is typically limited to two-dimensional (2D) analysis. In order to obtain a complete tissue deformation analysis, it is necessary to extend the implementation of the coupled filtering method to three-dimensional (3D) analysis. However, doing so will require billions of times more operations than the 2D implementation, which makes the 3D implementation impractical for real-world applications.

[0004] Therefore, a need exists to provide accelerated coupled filtering methods and systems for tissue deformation analysis.SUMMARY

[0005] According to a first aspect, there is provided an accelerated coupled filtering method for tissue deformation analysis, comprising: initializing, by a first processing device, a plurality of first motion parameters; applying, by a second processing device, a first filter on a pre-deformation image of a tissue to obtain a filtered pre-deformation image, wherein the second processing device is different from the first processing device; applying, by the second processing device, a second filter on a post-deformation image of the tissue to obtain a filtered post-deformation image, wherein the filtered pre-deformation image and the filtered post-deformation image are correlated by a first motion matrix comprising the plurality of first motion parameters, and wherein the plurality of first motion parameters comprise at least three fundamental first motion parameters, each of the three fundamental first motion parameters represents movement of the tissue along an axial direction relative to the axial direction, an elevational direction and a lateral direction, respectively; estimating, by the second processing device, a respective value for each of the plurality of first motion parameters, wherein each of the estimated respective values represents a difference between the pre-deformation and post-deformation filtered images; and in response to determining that at least one of the estimated respective values meets at least one predefined criterion, updating the at least one of the estimated respective values as an optimal value for at least one corresponding first motion parameter of the plurality of first motion parameters.

[0006] According to a second aspect, there is provided a system for performing an accelerated coupled filtering method for tissue deformation analysis, the system comprising a first processing device and a second processing device, wherein the second processing device is different from the first processing device, and wherein the first processing device is configured to: initialize a plurality of first motion parameters; and wherein the second processing device is configured to: apply a first filter on a pre-deformation image of a tissue to obtain a filtered pre-deformation image; apply a second filter on a post-deformation image of the tissue to obtain a filtered post-deformation image, wherein the filtered pre-deformation image and the filtered post-deformation image are correlated by a first motion matrix comprising the plurality of first motion parameters, and wherein the plurality of first motion parameters comprise at least three fundamental first motion parameters, each of the three fundamental first motion parameters represents movement of the tissue along an axial direction relative to the axial direction, an elevational direction and a lateral direction, respectively; estimate a respective value for each of the plurality of first motion parameters, wherein each of the estimated respective values represent a difference between the pre-deformation and post-deformation filtered images; and in response to determining that at least one of the estimated respective values meets at least one predefined criterion, update the at least one of the estimated respective values as an optimal value for at least one corresponding first motion parameter of the plurality of first motion parameters.BRIEF DESCRIPTION OF THE DRAWINGS

[0007] Embodiments of the invention will be better understood and readily apparent to one of ordinary skill in the art from the following written description, by way of example only, and in conjunction with the drawings, in which:

[0008] Fig. 1 shows a flowchart illustrating a conventional tissue deformation analysis process.

[0009] Fig. 2A (left) shows a visual illustration of a coarse-to-fine search strategy and Fig. 2A (right) shows a flowchart illustrating the coarse-to-fine search strategy of a conventional “baseline” coupled filtering method.

[0010] Fig. 2B shows a flowchart illustrating a detailed search process of a conventional “baseline” coupled filtering method.

[0011] Fig. 3 shows a table listing parameters used in FastCF, in accordance with an embodiment.

[0012] Fig. 4 shows a flowchart illustrating (a) a conventional process of converting a Radio Frequency (RF) image to B-mode image and (b) a process of a RF image to B-mode image in accordance with an embodiment.

[0013] Fig. 5 shows an example scenario of a post-processing step, in accordance with an embodiment.

[0014] Fig. 6 shows a flowchart illustrating an overview of an accelerated coupled filtering method ( “FastCF” ) , in accordance with an embodiment.

[0015] Fig. 7 shows a schematic diagram of an exemplary system for performing an accelerated coupled filtering method for tissue deformation analysis, in accordance with an embodiment.

[0016] Fig. 8 shows a flowchart illustrating an exemplary workflow of an accelerated coupled filtering method for tissue deformation analysis, in accordance with an embodiment.

[0017] Fig. 9 shows an exemplary model of an elastic tissue, in accordance with an embodiment.

[0018] Fig. 10 shows a table listing the time complexity, space complexity and estimated Floating-point Operations per second (FLOPs) for 2D and 3D implementations of various coupled filtering methods, in accordance with an embodiment.

[0019] Fig. 11 shows a table listing parameters used for the time complexity and space complexity analysis for the 2D and 3D implementations of the various coupled filtering methods, in accordance with an embodiment.

[0020] Fig. 12 shows another table listing detailed time complexity analysis for the 3D implementation of various coupled filtering methods, in accordance with an embodiment.

[0021] Fig. 13 shows a table listing the square root of mean square error (SRMSE) and computation time of FastCF performed on 3D simulated dataset, in accordance with an embodiment.

[0022] Fig. 14 shows color maps (in grayscale) illustrating results obtained from performing FastCF on the 3D simulated dataset, in accordance with an embodiment.

[0023] Fig. 15 shows a table listing computation time of different functions used in the FastCF process on the 3D simulated dataset, in accordance with an embodiment.

[0024] Fig. 16 shows color maps (in grayscale) illustrating results obtained from performing FastCF on 3D real data set, in accordance with an embodiment.

[0025] Fig. 17 shows a color map (in grayscale) of (a) baseline method and (b) FastCF performed on 2D simulated dataset, in accordance with an embodiment.

[0026] Fig. 18 shows graphs illustrating Normalized Correlation Coefficient (NCC) of images based on different motion parameters used in FastCF, in accordance with an embodiment.

[0027] Fig. 19 shows a table listing SRMSE of FastCF under different setting options, in accordance with an embodiment.

[0028] Fig. 20, comprising Fig. 20 (a) , Fig. 20 (b) and Fig. 20 (c) , shows histograms illustrating results of the accelerated coupled filtering method for tissue deformation analysis of Fig. 8.

[0029] Fig. 21 shows a schematic diagram of an example of a computing device used to realise a system for performing an accelerated coupled filtering method for tissue deformation analysis, in accordance with an embodiment.DETAILED DESCRIPTION

[0030] Embodiments of the present invention will be described, by way of example only, with reference to the drawings. Like reference numerals and characters in the drawings refer to like elements or equivalents.

[0031] Some portions of the description which follows are explicitly or implicitly presented in terms of algorithms and functional or symbolic representations of operations on data within a computer memory. These algorithmic descriptions and functional or symbolic representations are the means used by those skilled in the data processing arts to convey most effectively the substance of their work to others skilled in the art. An algorithm is here, and generally, conceived to be a self-consistent sequence of steps leading to a desired result. The steps are those requiring physical manipulations of physical quantities, such as electrical, magnetic or optical signals capable of being stored, transferred, combined, compared, and otherwise manipulated.

[0032] Unless specifically stated otherwise, and as apparent from the following, it will be appreciated that throughout the present specification, discussions utilizing terms such as “scanning” , “calculating” , “determining” , “replacing” , “generating” , “initializing” , “outputting” , or the like, refer to the action and processes of a computer system, or similar electronic device, that manipulates and transforms data represented as physical quantities within the computer system into other data similarly represented as physical quantities within the computer system or other information storage, transmission or display devices.

[0033] The present specification also discloses apparatus for performing the operations of the methods. Such apparatus may be specially constructed for the required purposes, or may comprise a computer or other device selectively activated or reconfigured by a computer program stored in the computer. The algorithms and displays presented herein are not inherently related to any particular computer or other apparatus. Various machines may be used with programs in accordance with the teachings herein. Alternatively, the construction of more specialized apparatus to perform the required method steps may be appropriate. The structure of a conventional computer will appear from the description below.

[0034] In addition, the present specification also implicitly discloses a computer program, in that it would be apparent to the person skilled in the art that the individual steps of the method described herein may be put into effect by computer code. The computer program is not intended to be limited to any particular programming language and implementation thereof. It will be appreciated that a variety of programming languages and coding thereof may be used to implement the teachings of the disclosure contained herein. Moreover, the computer program is not intended to be limited to any particular control flow. There are many other variants of the computer program, which can use different control flows without departing from the spirit or scope of the invention.

[0035] Furthermore, one or more of the steps of the computer program may be performed in parallel rather than sequentially. Such a computer program may be stored on any computer readable medium. The computer readable medium may include storage devices such as magnetic or optical disks, memory chips, or other storage devices suitable for interfacing with a computer. The computer readable medium may also include a hard-wired medium such as exemplified in the Internet system, or wireless medium such as exemplified in the GSM, GPRS, 3G or 4G mobile telephone systems, as well as other wireless systems such as Bluetooth, ZigBee, Wi-Fi. The computer program when loaded and executed on such a computer effectively results in an apparatus that implements the steps of the preferred method.

[0036] The present invention may also be implemented as hardware modules. More particularly, in the hardware sense, a module is a functional hardware unit designed for use with other components or modules. For example, a module may be implemented using discrete electronic components, or it can form a portion of an entire electronic circuit such as an Application Specific Integrated Circuit (ASIC) or Field Programmable Gate Array (FPGA) . Numerous other possibilities exist. Those skilled in the art will appreciate that the system can also be implemented as a combination of hardware and software modules.

[0037] In the following description, the term “module” can refer to software, a hardware element, or a combination of both.

[0038] An Application Programming Interface (API) enables software and applications to communicate with each other. It is a software-to-software interface that allows for separate parties to communicate with each other without any previous user knowledge or intervention. In general terms, it is a set of clearly defined methods of communication between various software components.

[0039] This specification uses the term “configured to” in connection with systems, apparatus, and computer program components. For a system of one or more computers to be configured to perform particular operations or actions means that the system has installed on its software, firmware, hardware, or a combination of them that in operation cause the system to perform the operations or actions. For one or more computer programs to be configured to perform particular operations or actions means that the one or more programs include instructions that, when executed by data processing apparatus, cause the apparatus to perform the operations or actions. For special-purpose logic circuitry to be configured to perform particular operations or actions means that the circuitry has electronic logic that performs the operations or actions.

[0040] As used herein, the term "processing device" refers to any hardware or system configured to perform computational tasks.

[0041] By way of examples, malignant tumours, cirrhotic liver, and dead heart tissues are significantly stiffer than healthy tissues. As such, stiffness can reveal pathologies, which led to doctors performing tissue deformation analysis to track the deformation induced by internal or external forces so as to deduce the stiffness distribution to facilitate diagnosis. Fig. 1 shows a flowchart illustrating a typical process of the tissue deformation analysis. For accurate tissue deformation analysis, coupled filtering method was proposed, which applies two different but coupled filters to images taken before and after tissue deformation, respectively, to compensate for both complex motion, i.e., feature motion decorrelation, and interference of echo waves. Adopting the coupled filtering method can result in a much higher accuracy than other relevant methods, especially for large tissue deformation. However, such method may require re-filtering of the images for different motions by exhaustively searching through all possible discrete motions. Consequently, such a heavy computational load makes only implementation on 2D images feasible in practical use, which may take around 45 hours to analyze a pair of images sized 201x1001 in MATLAB on the CPU and around 20 minutes on the FPGA.

[0042] Further, as tissues move in all directions, tissue deformation analysis is preferably performed on 3D images so that all motions can be examined to achieve high accuracy. Due to the exhaustive search for all possible affine motion, the 3D implementation of the coupled filtering method is estimated to perform billions more calculations than the 2D implementation. Therefore, the present disclosure seeks to improve the algorithmic complexity of the coupled filtering method. Specifically, the present disclosure provides a method that can accelerate the coupled filtering method while maintaining a similar accuracy (hereinafter referred to as “FastCF” ) . The main functions of FastCF are as follows: 1) approximating the motion can reduce the required number of times that the coupled filters are applied, 2) employing envelope detection can reduce the search space of the motion parameters by heuristically refining the search for increased speed and 3) adding a quick post-processing step can compensate possible inaccuracies introduced by the approximation. In the following section, the conventional coupled filtering method is briefly described, of which this method serves as a foundation to FastCF. Coupled Filtering Method and “Baseline Method”

[0043] The coupled filtering method can be, for example, designed to compensate for feature motion decorrelation. The ultrasound image before deformation, IA (X) , can be modeled as: IA (X) =Z (X) *H (X) ,                  (1) where Z (X) represents scatterers, * denotes convolution, and H (X) represents the point  spread function of the imaging system. Further, Z (X) and H (X) can be modelled as follows: Here, P= { (a1, p1) , (a2, p2) , ..., (aN, pN) }         (3) where a1, …, aN denotes amplitudes of the scatterers and p1, …, pN denotes positions of  the scatterers, X,  denote image coordinates and the positions of scatterers, respectively, N is the number of scatterers,  where x, y, and z represent the lateral, elevational and axial directions, respectively, and U0= [0 0 uz] Tare parameters for the point spread function.

[0044] Given the motion model: X′=M (X+T) ,           (4) with affine motion matrix and displacement vector T= [tx ty tz] T.

[0045] The following equation can be derived, which describes the coupled filtering method: IA (X) *H (MX) =IB [M (X+T) ] *H (X)         (5) where IB (X) is the image after deformation.

[0046] Specifically, the coupled filtering method combines the filtering step based on (6) with block matching to compensate for feature motion decorrelation and therefore can achieve high accuracy in large tissue deformation analysis. Figs. 2A (right) and 2B illustrate the workflow of the coupled filtering method. Note that the motion model (5) can allow the search for T by simply sampling across the inversely mapped image rather than repeating the entire coupled filtering procedure. This concept can help save computational load, which is adopted in the “baseline method” (described in the later section) .

[0047] Further, the coupled filtering method can adopt a coarse-to-fine search strategy to reduce the computational cost. As shown in Fig. 2A, this strategy first focuses on a small number of points and generates a rough low-resolution output, namely the first scale / iteration, and then gradually refines the result by increasing the resolution and using smaller search steps in the subsequent scales / iterations. Fig. 2B illustrates a detailed process flow of the search step. The iteration repeats until all points are processed. In the present disclosure, the coupled filtering method described above together with the coarse-to-fine search strategy adopted is referred to as the “baseline method” . FastCF

[0048] In this section, the present disclosure describes the algorithmic optimizations of FastCF. The post-processing step that helps maintain the accuracy of the estimated strain is also described. FastCF Parameters

[0049] Firstly, three search spaces are defined as follows: where denotes the coordinates of voxels in a  block.  and Mrange,  denote the range and step of the coordinates of points of interest, M, and T respectively;  denotes the element-wise less than operator;  is element-wise multiplication;  is the set of a vector sized 3-by-1;  is the set of a vector sized 3-by-3;  is the set of integers; I3 is the 3×3 identity matrix.

[0050] In further details,  denotes the points of interest. Theoretically, each point in the original image is only linked to another position in the deformed image (post-deformation image) . Hence, only two vectors are required to describe these two corresponding positions. Alternatively, the original position (i.e., V) and a displacement vector can also be used to describe these two corresponding positions. However, in order to determine if these two points matches, simply comparing the voxel value of these two points may not be enough because there may be other factors, such as noise or feature motion decorrelation, which may interfere with the accuracy of the result. As such, the voxel value of voxels surrounding the points of interest can be taken into consideration, i.e., matching a block of voxels in the original image (pre-deformation image) to a block of voxels in the deformed image. In this regard, a simple displacement vector may not be sufficient to account for such huge number of voxels (e.g., 200 voxels in each block) . As a result, a more complex motion model, i.e., the affine transformation, X′=M (X+T) , is adopted.

[0051] Next, the present disclosure describes how the search space can be reduced in order to speed up the search. From Equations (Error! Reference source not found. ) and (6) , the following Equation (7) can be obtained: IB [M (X+T) ] =Z (X) *H (MX)           (6)

[0052] Comparing the pre-deformation image as described in Equation (Error! Reference source not found. ) and the post-deformation image as described in Equation (6) , it can be seen that the key difference, among other things (if any) , lies in H (X) and H (MX) . Further, the point spread function H (X) as described in Equation (2) has a Gaussian shape in both x and y directions but is a Gaussian-weighted cosine function (e.g., Gabor function) in the z direction. Therefore, by comparing H (X) and H (MX) , it can be observed that changes in all three directions can cause interference between two Gaussian envelopes, but only changes in the z direction would cause significant interference between the peak and trough of the Gaussian-weighted cosine function, which can lead to larger decorrelation between the pre-deformation image (as described in Equation (Error! Reference source not found. ) ) and the post-deformation image (as described in Equation (6) ) and, in turn, degrades the accuracy of estimated stiffness.

[0053] However, if it is defined that the observation can be summarized as follows: H (McX) ≈H (MX)      (7) when |mxx-1|+|mxy|+|mxz|+|myx|+|myy-1|+|myz|≤0.6,

[0054] It may be worth reiterating that Mc and |mxx-1|+|mxy|+|mxz|+|myx|+|myy-1|+|myz|≤0.6 are formulated under the assumption that only changes in the z direction would cause significant interference between the peak and trough of the Gaussian-weighted cosine function. In other words, changes in the motion parameters, mxx, mxy, mxz, myx, myy and myz, may be negligible as long as they meet the above defined conditions. Additionally, under the above defined conditions, it should be noted that

[0055] In further details regarding Equation (8) , the following equations are first formulated based on both ends of Equation (6) : I1-ori (X) =IA (X) *H (MX)  (9) I2-ori (X) =IB [M (X+T) ] *H (X)   (10)

[0056] However, performing convolution for all possible Ms based on the above two Equations (9) and (10) can take too much time. Therefore, in order to reduce the amount of time taken to perform convolution, the following three Equations are further formulated: I1 (X) =IA (X) *H (McX)   (11) I2 (X) =I2t [M (X+T) ]   (13)

[0057] Using Equations (11) and (12) , I1 (X) and I2t (X) only need to be calculated when Mc changes. As Mc only has 3 independent variables, i.e., mzx, mzy and mzz, calculation time can be significantly reduced as compared to calculating M, which, at the maximum, has 8 independent variables (9 in total, minus one due to the assumption of tissue incompressibility) . Further, calculation of I2t [M (X+T) ] , i.e., Equation (13) , can be performed by re-sampling the image. As a result, the time complexity of calculating Equation (13) is whereas the time complexity of convolution is Therefore, it is evident that adopting Equations (11) to (13) is more efficient.

[0058] Further, based on Equations (12) and (13) , the following Equation (14) can be derived:

[0059] As described above,  when under the defined conditions in relation with Equation (8) . Therefore, the following relationship can be obtained:

[0060] Similarly, based on Equation (8) , the following Equation (16) can be derived: I1 (X) =IA (X) *H (McX) ≈IA (X) *H (MX) =I1-ori (X)   (16) where S is defined as the volume of the imaging area. This condition defines that the  number of scatterers used in the simulation should reach around 10 per resolution cell to form the speckle in ultrasound images.

[0061] As Equation (6) defines that IA (X) *H (MX) =IB [M (X+T) ] *H (X) , then I1 (X) ≈I2 (X) .

[0062] In summary, based on Equations (14) to (16) , the following relationship can be derived:

[0063] In particular and as would be appreciated by the person skilled in the art, the approximate Equation (17) indicates that after the calculation of IA (X) *H (McX) and  if the changes in the other six motion parameters fulfil the above defined condition, i.e., less than the threshold defined above, there is no need to recompute the convolution. Advantageously, as convolution is more time-consuming than computing I2t [M (X+T) ] , this implementation can help save time. Further, the filtering scheme in different stages of coarse-to-fine iterations can be adjusted. For example, during the last two iterations, or namely the last two scales, the filtering step would account for over 90%of the total computation load. Thus, the user may apply the filtering step only once without compromising on the accuracy of the results. Specifically, since the searching step Mstep is designed to be small in comparison with the search space in these scales, such a small shift would not degrade the result significantly.

[0064] Also, in example embodiments, three motion parameters, for example mxy, myy, and mzy, can be omitted from the motion model as described in Equation (5) . This omission can be done by flattening the 3D block along the y (elevational) direction, which makes these three variables / motion parameters meaningless when calculating M (X+T) . Originally, the incompressibility constraint already omits one motion parameter from the total number of motion parameters to be searched, i.e., reduced from a total of 9 motion parameters to 8 motion parameters. However, in the above example embodiments, three motion parameters are eliminated, i.e., mxy, myy, and mzy, which makes the incompressibility constraint to no longer take effect. Hence, mxy, myy, and mzy are omitted from a total number of 9 motion parameters (instead of 8) . Consequently, this implementation can advantageously accelerate the computation by reducing the total motion parameters to be searched by two motion parameters from M, i.e., from 8 motion parameters (incompressibility constraint is in effect) to 6 motion parameters. It would also be appreciated by the person skilled in the art that more motion parameters may be removed, depending on the desired accuracy. Additionally, in the above example embodiments, the axial direction (z direction) is retained because it has the highest resolution. Further, the elevational direction (y direction) is omitted / flattened because it usually has lower resolution than the lateral direction. In other words, retaining the axial and lateral directions may ensure accurate results to be obtained.

[0065] Putting it all together, Mc can be defined as the matrix used in the filtering step as follows:

[0066] Algorithm 1 below presents the algorithm with the updated filtering step. The envelope function will be described in the next section. In Algorithm 1, multiple loop sequences are defined for different search scales. Specifically, with reference to Fig. 3, in the first scale / iteration (scale / iteration =1) , the search range of T, i.e.,  is relatively large, whereas the search space of M, i.e.,  is the same for all the points of interest, so it is worth calculating on the entire image. On the other hand, the subsequent scales / iterations (scale / iteration > 1) have smaller and different points have different  so it is suitable for calculating only part of the image. Envelope Detection

[0067] After the filtering step as described above, the images may retain wavy shapes, and their envelopes can contain useful information for deformation analysis. Therefore, embodiments of the present disclosure provide a method to detect these envelopes and aim to take advantage of the image envelope (s) in the filtering step such that the search can proceed at much larger steps at the coarsest scale. In more details, Radio-Frequency (RF) signal is the raw signal coming out of the ultrasound sensor, and B-mode signal is easier to interpretate. Embodiments of the present disclosure provide a different workflow to detect the envelope (s) of filtered RF signal and convert it into filtered B-mode signal (s) , as shown in Fig. 4, which allows the search steps to be enlarged without causing severe peak hopping issue in RF signal-based motion tracking. Specifically, with reference to Fig. 4, a filtering step can be applied to an RF image, e.g., pre-deformation image or post-deformation image. Thereafter, detection of envelope (s) can be performed to the filtered RF image. Finally, a spatial transformation step can be performed to obtain a Filtered B-mode image. The order of performing the above-mentioned process steps can be changed. Advantageously, performing the spatial transformation step after filtering the RF image can improve computation speed.

[0068] Additionally, the original RF signal can have a very large spatial frequency in the axial direction, e.g. for a 3MHz ultrasound signal and a sampling speed of 15.4 MHz, peak and trough can be found in every 5.1 voxels. This can mean that, with a slight shift in the M and T, the resulting image I (MX+T) may look different, and it can be easy to mismatch the peaks and troughs. Hence, the search steps of M and T are set to be relatively small so that the blocks can be correctly matched by combining the pattens of an entire block. With envelope detection, the spatial frequency of the wavy signal can be decreased by around 5 to 7 times as shown in Fig. 4, e.g., one peak and one trough for every 30 voxels. This means that the original search step of M and T can be increased by around 6 times, and it would still enable the correct matching of peaks and troughs while mitigating the peak hopping issue. Consequently, a rough search can be applied at the beginning and the result can be finalized with smaller search steps thereafter. Importantly, it would be appreciated by the person skilled in the art that the above conversion process can be applied together with the coarse-to-fine multi-scale strategy for further acceleration. The envelope is the L1-norm of the Hilbert transformed signal. The algorithm of the updated coarse-to-fine strategy is shown in Algorithm 2.

[0069] It may be worth mentioning that the flattening is performed on the blocks (or one or more voxels) of the images, rather than on the images. In other words, when using a similarity metric to measure the similarity between two corresponding points in pre-and post-deformation images, only a flattened region of interest is considered, i.e., blocks including X*1*Z number of voxels, in both the pre-and post-deformation images. Specifically, when calculating the similarity metric, only the voxels in the 2D blocks are considered. However, in order to accurately calculate the values of these 2D blocks, it may still be necessary to consider the voxels surrounding the voxels of these 2D blocks as the voxels in other elevational planes may affect the accuracy of the result. Post-Processing

[0070] Peak hopping is common in RF image-based tissue deformation analysis, and it can degrade the motion tracking results. To address this issue, embodiments of the present disclosure provide a post-processing step based on majority voting. Particularly, it is observed that the axial displacement, i.e., Dz, is more accurate than the other two directions. Consequently, based on this observation, limits can be imposed on the gradients of Dz in three directions. For example, by defining two adjacent points P,  with P-Q= [Vstep, x 0 0] T, and the displacement Dopt (V) =Mopt (V) (V+Topt (V) ) -V, with the Lagrangian motion model, the following limits can be derived: where [mzx, min, mzx, max] is the search range of mzx.

[0071] Similar limits can also be established along the other two directions. With reference to Fig. 5, any gradient exceeding these theoretical limits can be recognized to be wrong, and the voxels surrounded by wrong gradients can be defined as wrong voxels. Subsequently, the M and T of the wrong voxels can be interpolated based on the correct voxels and let the finer scales refine the results. Usually less than 5%are wrong voxels, and other voxels are not affected. Also, the post-processing step can be designed to be light-weight and does not substantially slow down the overall computing time. The results of the post-processing step are discussed in the later section of this disclosure. The post-processing step algorithm is shown in Algorithm 3 below. FastCF Algorithm

[0072] The overall algorithm of FastCF is presented in Algorithm 4 and described in this section. When searching for the optimal T in Algorithm 1, the sum-tables method can be adopted. In example embodiments, the CPU implementation is coded in MATLAB. In the first scale / iteration, the computation is parallelized after calculating I1 (X) and I2temp (X) ; in the remaining scales / iterations, the computation of all Vs is parallelized. Further, the MATLAB code for the interpolation and searching for the best T with the C++ MEX functions is also re-implemented. Pointer and dimension permutation are also adopted. In the post-processing step, the graph can be stored in the form of 3D matrix, of which elements with odd coordinates represent the vertices as shown in Fig. 5. Each vertex represents one voxel in  and can be connected to its adjacent vertices only when the difference of their Dopt falls into the theoretical limit. The MATLAB built-in function bwlabeln () is then utilized, which helps find connected components in binary image, to accomplish the relevant graph operations.

[0073] Figs. 6 provide an overview of an exemplary process flow of the method, i.e., “FastCF” , described in the present disclosure. In other words, embodiments described in the present disclosure provide an accelerated coupled filtering method for tissue deformation analysis. Graphical Processing Unit (GPU)

[0074] The following paragraphs provide a brief background of GPU and programming using the GPU. GPU is an acceleration hardware that is designed for computational-intensive tasks. Although primarily used to process and render graphics, the GPU’s ability to process huge amounts of data and perform parallel computations make it useful for fields such as Artificial Intelligence, image processing, computing acceleration, etc. A typical GPU can have thousands of computing cores. Every few computing cores may be grouped as a warp and can work together by communicating with one another. However, if the computing cores in a warp need to execute different operations, this can create a warp divergence and the warp can only execute one branch at a time, thereby degrading the throughput. In embodiments described in the present disclosure, workflow of the algorithm has been designed to ensure that the thousands of computing cores can work concurrently while creating few warp divergence. As a non-limiting example, on Nvidia GPU GeForce RTX 2080 Ti, a warp size can be 32 and a total number of computing cores can be 4352. A person skilled in the art would appreciate that GPUs from other brands can be used.

[0075] In a GPU program, the computing cores may be organized by grids and blocks. A layer -thread block clusters may also be present in some GPUs. Implementations of the invention may only require two levels, hence the thread block clusters may not be utilized. A person skilled in the art would appreciate that in some embodiments, the thread block clusters may be utilized. Each grid of computing cores may contain several blocks, and the blocks can be organized into 1D, 2D or 3D arrays. Each block may contain several threads, and the threads can be organized into 1D, 2D or 3D arrays. The programming model may assume that the GPU is a co-processor of a CPU. To run a program on the GPU, kernel functions designed for multi-threading on thousands of computing cores may be executed, and defined sizes of grids and blocks may be passed as a variable for the GPU to schedule tasks.

[0076] Both the CPU and the GPU have their own memories to store data during runtime. In the GPU program, the former is termed as host memory and the latter is termed as device memory. As the host memory and the device memory are two separate memories, communication between the memories may go through PCIe Bus, which is slower than both memories. Therefore, in embodiments of the present invention, most of the computing tasks may be implemented on the GPU so that the impact of limited bandwidth between the CPU and the GPU can be reduced, except for initializing parameters that may be implemented on the CPU. Initializing parameters is one of the few tasks that may be implemented on the CPU because it may include many branches and small loops and may create excessive warp divergence. Moreover, data size may be too small to take advantage of thousands of GPU cores. Advantageously, using the GPU for the remaining computational-intensive tasks can improve processing speed. For example, a speed up factor of 10 is achieved and results of tissue deformation analysis may be obtained in 2.5 minutes. Further, the GPU with equivalent or higher performance as compared to the CPU can cost less than the CPU. This lowers the barriers for a clinician to adopt the method and system.

[0077] There may be long latency to access device memory. Retrieving one byte of data from a memory can take hundreds of computing cycles, which may lead to time degradation and inefficiency. To circumvent the long latency, rules may be implemented on embodiments of the invention. The rules may include but not limited to 1) principle of locality; 2) accessing data in a continuous manner so that a memory access operation in a warp can be coalesced; and 3) allocate more concurrently running jobs than the actual number of computing cores, so that when there exists extra registers, the scheduler can run several threads at the same time and perform a relatively fast context switching to hide the latency of memory access.

[0078] A workflow of the GPU-accelerated implementation according to an embodiment of the invention, is presented in algorithm 5. The computation, unless specifically mentioned, can be performed by several GPU kernel functions with more details provided below.

[0079] The Fourier transform can be the most time-consuming task in the accelerated coupled filtering method. As a non-limiting example, CUDA built-in library cuFFT (Nvidia CUDA Fast Fourier Transform) can be used. For the first scale, since Fourier transform only needs to be applied a few times to a large image, the library may be called directly. For the following scales, the Fourier transform may be applied more times on smaller image blocks. cuFFT batch plan that allows application of the Fourier transform on a batch of images at the same time can be used.

[0080] In implementations of the invention, the size of the FFT may be equal to the size of the block plus the size of the point thread function minus one. The size of the point spread function may be the one when its envelope equals to 1 / 100 of its amplitude. The cuFFT can achieve higher efficiency when dealing with images sizes that can be represented in a form of 2a×3b×5c×7d. The smaller the prime factor, the faster the speed that the cuFFT can achieve. In other words, the image size that can be represented by 2a can be processed at the fastest speed. Therefore, the image size can be adjusted by padding zeros for faster processing.

[0081] Calculating sum table and correlation coefficient may be a second most computational intensive task. Aside from using the sum table in the accelerated coupled filtering method to prevent redundant computation, two means to improve speed of computation in the GPU implementation can be used and algorithm 6 is provided below to illustrate the workflow.

[0082] First, the sum table can be calculated block by block. For example, for the first scale, there are many voxels that need to compute the sum table. To make better use of L1 cache, the computation can be separated into several tiles.

[0083] Second, dimensions can be permuted after calculating the sum table. The image may be originally stored as a 3D array with dimensions (x, y, z) . However, when the block is flattened in the accelerated coupled filtering method, it is flattened in the y direction. As such, the sum table may only need to be built along the (x, z) directions. Therefore, when storing the sum table, the dimensions to (y, x, z) may be permuted so that one sum table can be stored continuously along linear addresses of the device memory. The principle of locality may be followed, and the program can run faster.

[0084] Moving on to the size of grids and blocks, for the first scale, a number of blocks (i.e. windows) can be directly passed at different locations as the grid size. For the remaining scales, since the number of blocks grows exponentially and can cause memory overflow, the grid size may be set as 1024, which is the batch size that the program can process at the same time and can be adjusted to save memory usage. The block size can range from 32～1024 depending on the computational load.

[0085] For remaining tasks such as multiplying a given number to normalize the result, they may not need to be done in a block wise manner. For example, the grid and block size can be set as 68 and 256 respectively, whose multiplication equals 4 times the number of computing cores on GPU 2080Ti. This can be slightly faster than smaller grid sizes as more blocks may allow for quick context switching to hide the memory access latency.

[0086] Further, a few tricks may be introduced. For example, when finding the largest NCC, the current best result may be stored in the shared memory and a reduction tree may be adopted for computing the final result. The allocated space for Fourier transformation of images can be reused to detect envelopes to save memory usage. Further, double buffering may be adopted to facilitate the transition between different scales.

[0087] The accelerated coupled filtering method for tissue deformation analysis can be implemented with system 700 shown in Fig. 7, which shows a schematic diagram of the system 700 for performing the accelerated coupled filtering method for tissue deformation analysis. The system 700 can include one or more processing devices 702. In embodiments of the present disclosure, the system 700 may be communicatively coupled to an ultrasound system 704 to receive and send information / data, e.g., ultrasound images, etc. However, the person skilled in the art would appreciate that the system 700 can be communicatively coupled to other suitable system or device for performing the method, depending on the application.

[0088] With reference to Fig. 8, the method 800 can include the following steps:

[0089] Step 802: initializing, by a first processing device, a plurality of first motion parameters.

[0090] Step 804: applying, by a second processing device, a first filter on a pre-deformation image of a tissue to obtain a filtered pre-deformation image, wherein the second processing device is different from the first processing device.

[0091] Step 806: applying, by the second processing device, a second filter on a post-deformation image of the tissue to obtain a filtered post-deformation image.

[0092] The filtered pre-deformation image and the filtered post-deformation image are correlated by a first motion matrix including the plurality of first motion parameters. The plurality of first motion parameters include at least three fundamental first motion parameters. Each of the three fundamental first motion parameters represents movement of the tissue along an axial direction relative to the axial direction, an elevational direction and a lateral direction, respectively.

[0093] Step 808: estimating, by the second processing device, a respective value for each of the plurality of first motion parameters. Each of the estimated respective values represents a difference between the pre-deformation and post-deformation filtered images.

[0094] Step 810: in response to determining that at least one of the estimated respective values meets at least one predefined criterion, updating the at least one of the estimated respective values as an optimal value for at least one corresponding first motion parameter of the plurality of first motion parameters.

[0095] In some implementations, the initializing step of step 802 may include loading two images, calculating all required parameters, allocating and freeing spaces on both the first processing device and the second processing device.

[0096] It would be appreciated by the person skilled in the art that, based on the term “corresponding” , it is defined, by way of example, that if the estimated respective value of only the fundamental motion parameter representing movement of the tissue along the axial direction relative to the lateral direction meets the at least one predefined criterion, the respective value of only such fundamental motion parameter will be updated. The remaining fundamental motion parameters will not be updated.

[0097] In embodiments of the present disclosure, the first processing device may be a Central Processing Unit (CPU) and the second processing device may be a Graphical Processing Unit (GPU) .

[0098] The method 800 can be used for three-dimension radio-frequency images. The method 800 can also be adjusted to be used for two-dimensional and common B-mode ultrasound images.

[0099] In embodiments of the present disclosure, the method 800 can further include executing, by the GPU, at least one kernel function associated with multi-threading on a plurality of computing cores. As would be appreciated, the at least one kernel function can be dozens of kernel functions. As a non-limiting example, the at least one kernel function can be 35 kernel functions. In some implementations, the method 800 can further include executing, by the GPU, common techniques. As non-limiting examples, the common techniques may comprise reduction tree, tiling, data permutation, double buffering and reusing allocated space. The common techniques may be adopted for faster execution speed and lower memory usage on the GPU.

[0100] In embodiments of the present disclosure, the method 800 can include a step of detecting, by the second processing device, one or more envelopes present in the filtered pre-deformation image and filtered post-deformation image to obtain a filtered B-mode pre-deformation image and a filtered B-mode post-deformation image.

[0101] In embodiments of the present disclosure, the applying steps 802 and 804 can include convolving a first point spread function and second point spread function of an ultrasound system with the pre-deformation image and the post-deformation image, respectively. The first point spread function can be a modified version of the second point spread function.

[0102] Additionally, after applying step 804, the method 800 can include a step of spatially transforming, by the second processing device, the post-deformation image based at least on a second motion matrix. The first motion matrix can be a modified version of the second motion matrix.

[0103] In embodiments of the present disclosure, the filtered pre-deformation image and filtered post-deformation image can be 3-dimensional (3D) images. The person skilled in the art would be able to appreciate that the pre-deformation, post-deformation, filtered pre-deformation and filtered post-deformation images are not solely limited to being 3D images and can encompass other image format such as 2-dimensional images, etc., depending on a user’s preference and / or the application. Further, estimating the respective value can include the following sub-step.

[0104] Sub-step A: flattening, by the second processing device, a first and second plurality of voxels surrounding each of one or more pre-deformation points of interest in the filtered pre-deformation image and each of one or more post-deformation points of interest in the filtered post-deformation image, respectively. The flattening can be performed along the elevational direction of each of the first and second plurality of voxels, which can reduce the plurality of first and second motion parameters.

[0105] In embodiments of the present disclosure, the estimating step 806 can be performed in a plurality of iterations. In a first iteration of the plurality of iterations, the plurality of first motion parameters can include two of the three fundamental first motion parameters representing movement of the tissue along the axial direction relative to the axial direction and the lateral direction, respectively, and the remaining fundamental first motion parameter may be removed by the flattening. Between a second iteration to a final iteration of the plurality of iterations, the plurality of first motion parameters can include the two of the three fundamental first motion parameters, two lateral first motion parameters representing movement of the tissue along the lateral direction relative to the lateral direction and the axial direction, respectively, and two elevational first motion parameters representing movement of the tissue along the elevational direction relative to the lateral direction and the axial direction, respectively.

[0106] In embodiments of the present disclosure, the method 800 can further include executing, by the second processing device, one or more Fast Fourier Transform batch plans. The one or more Fast Fourier Transform batch plans may be configured to perform Fourier Transform on a plurality of images in the first iteration. Alternatively, the one or more Fast Fourier Transform batch plans may be configured to perform Fourier Transform on the first and second plurality of voxels surrounding each of the one or more pre-deformation points of interest in the filtered pre-deformation image and each of the one or more post-deformation points of interest in the filtered post-deformation image, concurrently. Size of the first and second plurality of voxels can be adjusted by padding zeros for faster executing efficiency on the second processing device.

[0107] A Fourier Transform kernel size can be adjusted so that the cuFFT library can use the hardware more efficiently while creating almost the same result.

[0108] In embodiments of the present disclosure, the second motion matrix can include a plurality of second motion parameters and the estimating step 806 can further include the following sub-steps.

[0109] Sub-step B: searching, by the second processing device, a value of each of the plurality of first motion parameters for each of the one or more pre-deformation points of interest in the filtered pre-deformation image.

[0110] Sub-step C: searching, by the second processing device, a value of each of the plurality of second motion parameters for each of the one or more post-deformation points of interest in the filtered post-deformation image.

[0111] Sub-step D: generating, by the second processing device, one or more pre-deformation voxels based on the searched value of each of the plurality of first motion parameters.

[0112] Sub-step E: generating, by the second processing device, one or more post-deformation voxels based on the searched value of each of the plurality of second motion parameters.

[0113] Sub-step F: calculating, by the second processing device, a similarity metric based on one or more first voxels surrounding the one or more pre-deformation voxels and one or more second voxels surrounding the one or more post-deformation voxels. The similarity metric can define a similarity between a target pre-deformation point of interest and a post-deformation point of interest corresponding to the target pre-deformation point of interest.

[0114] Sub-step G: determining, by the second processing device, if the target pre-deformation point of interest matches the corresponding post-deformation point of interest based on the calculated similarity metric.

[0115] Additionally, in embodiments of the present disclosure, the plurality of second motion parameters can include at least two fundamental second motion parameters, two lateral second motion parameters and two elevational motion parameters. The estimating step 806 can further include, in the first iteration, determining, by the second processing device, a value of each of the two of the three fundamental first motion parameters based on steps (A) to (G) . Between the second iteration to a third final iteration, determining, by the second processing device, an updated value of each of the two of the three fundamental first motion parameters based on steps (A) to (G) , and assigning, by the second processing device, a value of each of the two lateral second motion parameters and the two elevational second motion parameters searched in step (C) in a previous iteration as a respective value of each of the two lateral first motion parameters and the two elevational first motion parameters. In a final two iterations, assigning, by the second processing device, the determined updated value of each of the two of the three fundamental first motion parameters in the third final iteration as an optimal value of each of the two of the three fundamental first motion parameters, and assigning, by the second processing device, an updated value of each of the two lateral second motion parameters and the two elevational second motion parameters searched in step (C) in the third final iteration as an optimal value of each of the two lateral first motion parameters and the two elevational first motion parameters.

[0116] In embodiments of the present disclosure, the similarity metric may include any one of the following: Normalized Correlation Coefficient, Sum of Absolute Differences and Sum of Squared Differences. However, the person skilled in the art would readily appreciate that other suitable metrics may be used, depending on a user’s preference and / or application.

[0117] Furthermore, the method 800 can include the following steps: imposing, by the second processing device, an upper limit and a lower limit on a gradient of a displacement range of the fundamental first motion parameters along the axial direction. Values falling outside the displacement range can be determined to be a wrong value of the fundamental first motion parameters. In response to determining that one or more voxels include the wrong value of the fundamental first motion parameters, determining, by the second processing device, a corrected value of the fundamental first motion parameters of the one or more voxels. The corrected value can be determined by interpolating the wrong value based on a predetermined correct value of the fundamental first motion parameters of voxels surrounding the one or more voxels. As one would appreciate, the first processing device can be used to impose the upper limit and the lower limit on the gradient of the displacement range of the fundamental first motion parameters along the axial direction. Also, in response to determining that one or more voxels include the wrong value of the fundamental first motion parameters, the first processing device can be used to determine the corrected value of the fundamental first motion parameters of the one or more voxels. Evaluation Experimental Setup

[0118] In this section, the time performance and tracking accuracy of FastCF are evaluated. The experimental setup includes hardware, software, two simulated and one real dataset, a set of parameters, and an evaluation metric. The machine used in the experiment includes a workstation with CPU AMD Ryzen Threadripper PRO 3995WX (64 cores) with 512GB memory, a GPU server with 2 CPUs Intel Xeon Gold 6130 (2x16Cores) with 256 GB memory and 8 Nvidia RTX 2080Ti with 11 GB device memory each, and another GPU virtual machine with 8 cores from Intel Xeon Gold 6246R, 32GB memory and 1 Nvidia A100 40 GB PCIe. The CPU program utilizes less than 64 GB of memory with 64 MATLAB workers. The CPU code was implemented and tested in MATLAB R2021b on Windows server 2022, and the GPU code was developed in CUDA version 12.1 with common compilation options in Ubuntu 18.04. Dataset

[0119] To simulate a dataset with known ground truth, SolidWorks is used to model an elastic tissue. As shown in Fig. 9, the elastic tissue has a hard spherical area merged in a soft cuboid material. The model is then uniformly compressed by 5%in the axial direction and SolidWorks is used to simulate the movement under compression. Further, one pair of 3D ultrasound images is generated using simulation software Field II based on the simulation result and estimate the parameters of the point spread function H (X) . The RF format ultrasound image is stored, and the image size is 101x101x1001 voxels. Additionally, a real 3D dataset is collected with the Philips iE33 xMATRIX ultrasound system. A X3-1 transducer is adopted to acquire 3D B-mode images of a participant’s left ventricle, and the image size is 302×158×302 voxels. Two volumes are then manually selected during the contraction for illustration. Due to lack of access to the RF signal, an approximate version of the coupled filtering method is adopted by removing the filtering step as follows: IA (X) ≈IB [M (X+T) ]          (20)

[0120] Further, FastCF is also tested using 2D simulation dataset and the results are used for comparison with the baseline method. Parameter Setting

[0121] The parameters and values of the simulated 3D dataset are defined in Fig. 10. Further, Fig. 10 also shows the time complexity, space complexity and estimated Floating Points per Seconds (FLOPs) of the 2D and 3D datasets. In Fig. 10, (a) represents coupled filtering with exhaustive search, (b) represents coarse-to-fine search strategy, (c) represents approximation of matrix M, (d) represents envelope detection, and (e) represents sum table. The other parameters are shown in Fig. 3. Additionally, Fig. 11 shows the parameters used for the time complexity and space complexity analysis for the 2D and 3D implementations of the various coupled filtering methods. In the experiment, the search range of mxx and mzz is between -10%and 10%, and the search range of mxz, myx, myz, mzx is between -5%and 5%. These search ranges are sufficiently large for the dataset. Evaluation Metrics

[0122] In the experiment, the square root of the mean square error (SRMSE) of the axial strain is used to measure the accuracy of estimation result: where ε (W) =mzz-1 is the estimated axial strain for a certain block B, ε0 (B) denotes  the corresponding ground truth, and NB is the number of blocks.

[0123] In the MATLAB program, a timer “tic toc” is used to measure the running time and “mpiprofile” is used to measure the running time of different functions. Only the time initializing the MATLAB parallel pool is excluded. FLOPs Estimation and Complexity Analysis

[0124] Fig. 12 shows the detailed time complexity of all related methods. The present disclosure also demonstrates how the time complexity and FLOPs of the exhaustive search are calculated as an example.

[0125] For the filtering step, kernels sized are convolved with images sized which can be accelerated by the circular convolution property of the Fourier Transform. Specifically, Fast Fourier Transformation (FFT) is applied on the padded image and the result is multiplied with the transformed kernel, then inverse FFT is applied and the final result is obtained, i.e., 4 forward FFT and 2 inverse FFT to get two convolution results IA (X) *H (MX) and IB [M (X+T) ] *H (X) . As IA (X) and H (X) remain constant, they only need 2 forward FFT for the entire program. Meanwhile, 1D FFT on data sized 1xn takes 5nlog2 n times operations. 3D FFT shares similar definitions except three nested summations, and each summation has the same complexity as 1D FFT, thus 3D FFT takes 15n3 log2 n3 operations. Inverse FFT also has the same operations as FFT, and multiplication takes n3 operations. Thus, the total FLOPs of the filtering step are and the time complexity is

[0126] Further, block-matching algorithm is adopted to search for the best results. For the interpolation step, interpolation is only performed on image for each M, and T can be achieved by sliding the window in the interpolated image. This implementation shows 30 operations for each voxel, so the FLOPs in total are and time complexity is  As for metric calculation, for voxels in the output, search is performed through Ms and Ts for its surrounding block sized and 3 additions and 3 multiplications are performed on each voxel. Each metric also has few operations but is negligible compared to others. Therefore, the rough FLOPs of block matching are  and the time complexity is In all, the total time complexity is

[0127] Further, by adopting the approximation of M (as described above) and with reference to Fig. 12, the convolution can be reduced from times to and the time to calculate the affine motion I2t [M (X+T) ] can be reduced from to Consequently, the time complexity can drop to If the na blocks are simultaneously processed, the space complexity can become Furthermore, if envelope detection is applied, approximation of matrix M can reduce the number of times to calculate the affine motion from to a constant, so the main computation load shifts back to convolution and the time complexity drops to Results on 3D Simulated Dataset

[0128] On the simulated 3D dataset, the estimated FLOPs and complexities of different methods are shown in Fig. 10. Both methods of M approximation and envelope detection effectively reduce the load and the usage of the sum table method also helps. In all, FastCF is estimated to reduce the FLOPs by 4.21 billion times compared to the baseline method. The mean SRMSE and the computing time at each scale are shown in Fig. 13. The results are shown in Fig. 14. As shown in Fig. 13, the SRMSE constantly decreases as the search steps of M and T get smaller. Further, with reference to Fig. 14, the final result clearly shows the hard sphere in the section view (Figs. 14 (a) and 14 (c) ) , fulfilling the goal of tissue deformation analysis. When taking a close look at the running time of different scales, it is observed that, except for the first scale (which initializes the search) , each of the other scales exhibit an increase of the number of blocks by around eight times. As such, the gradual decrease of the search space is not enough to compensate for the significant increase of number of blocks, which led to the overall computation time increasing. Additionally, as described above, the filtering strategy at the fifth scale is adjusted. As such, the fifth scale takes shorter time than others. The total time is slightly longer than the summation of all six scales’ time because the interpolation between scales takes a few seconds.

[0129] Further, the computation time of different functions is analyzed, as shown in Fig. 15. In this analysis, the program is split into four main components, and the first two main components are further broken down. As shown in Fig. 15, the coupled filtering component is the most time-consuming, where the Fourier Transform relies on built-in libraries and cannot be externally optimized. Data preparation and frequency domain multiplication are likely limited by the memory speed. Further, the post-processing component takes a negligible amount of time as compared to the search process. The long computation time of other functions in the CPU implementation is mostly spent by MATLAB indexing. Results on 3D Real Dataset

[0130] For the real 3D dataset, the white area near the ventricle wall is manually segmented in advance. The result is shown in Fig. 16. The analysis takes FastCF around 295 seconds to generate a strain image sized 97×65×129 on the CPU. With reference to Fig. 16, the overall compression trend in the slices is clearly visible, and further interpretation of these results would require clinical information. This dataset was generated with a frame rate of only 8 Hz. In other words, embodiments of the present disclosure can generate decent results with the ultrasound images sampled at a lower framerate, e.g., 8 Hz. In contrast, traditional deformation analysis methods require a frame rate of at least 40Hz to 80Hz to generate decent results. Comparison with Baseline Method

[0131] For comparison purposes, FastCF and the baseline method are performed on 2D simulated data that contains a hard cylinder. This is because performing the baseline method on 3D data will take too much time. The ground truth is similar to Fig. 14 (a) . As shown in Fig. 17, the baseline method and FastCF can both identify the hard inclusions in the 2D dataset. Further, FastCF achieves an SRMSE of 1.11%whereas the baseline achieves an SRMSE of 1.40%. FastCF also finishes the analysis within around 12.9 seconds on the CPU, but the baseline method takes around 1.72 days on the same CPU and 20 minutes on FPGA. Verification of the Approximation of Matrix M

[0132] In this subsection, the effect of the three parameters in M (as defined in Equation (5) ) is studied. Firstly, the image model in Equation (1) is used and the scatterer-density requirement as defined in Equation (8) is applied to generate the pre and post deformation images with different motion matrices M. Then, 27 blocks are evenly selected and their normalized correlation coefficient (NCC) is averaged out to plot the curves between the NCC and the nine motion parameters in M separately. The result is shown on Fig. 18. With reference to Fig. 18, it is observed that the results are consistent with the Equation (8) , i.e., if there are only small shifts along mx* and my*, while mz* does not change, IA (X) and IB [M (X+T) ] still can achieve high correlation coefficient. Study of Search Space and Post-Processing

[0133] Further, in this experiment, four options were tested without post-processing, i.e. “veryslow” to “medium” , to study the effect of different search spaces of M and T. The running time of these options are 434.607s, 83.928s, 15.796s, 2.571s and 2.589s respectively. For comparison purposes, the FastCF, i.e. the “medium” setting with post-processing, is also included in this experiment. Fig. 19 shows the errors across six scales using five options. Compared to the “medium” setting without post-processing, the slower options achieve lower SRMSE but show diminishing returns. On the other hand, FastCF is able to achieve similar accuracy as other slower options and runs nearly as fast as the “medium” option. Implementing with GPU

[0134] The experiment shows that given the same input, the processing speed improved by 10 times and the cost of the system is at least 2 times cheaper, and newer model could provide faster processing speed while keeping the cost the same. In implementations with only the CPU with MATLAB and MEX functions, it takes 2589 seconds to complete the computation. On the other hand, in implementations using a GPU, with a program using only one RTX 2080Ti, the program takes 274 seconds to complete computation, i.e. 1 / 10 of the original time. In addition, the system cost with the GPU is significantly lower than that with only the CPU.

[0135] Testing is also done using a relatively more powerful GPU A100 40GB PCIe, and computation is completed within 135 seconds. As Nvidia RTX 2080Ti is launched in 2018, the test on A100 shows possible running time that can be achieved on newer GPUs, and it would still cost less compared to the workstation CPU. Further, when running 8 separate tasks on 8 different GPUs simultaneously on the GPU server, throughput is nearly enlarged by 8 times, which can be helpful in clinical environment.

[0136] Fig. 20, comprising Fig. 20 (a) , Fig. 20 (b) and Fig. 20 (c) , shows histograms illustrating results of the accelerated coupled filtering method for tissue deformation analysis of Fig. 8. Specifically, Fig. 20 shows that results from a CPU implementation and a GPU implementation of the accelerated coupled filtering method are almost the same. Fig. 20 (a) shows a histogram of CPU implementation results, specifically, estimated axial strain from the CPU implementation and Fig. 20 (b) shows a histogram of GPU implementation results on 2080Ti. Fig. 20 (c) shows a histogram of a difference between the CPU implementation results and the GPU implementation results, specifically, pointwise differences of the estimated axial strain. As can be seen, the majority of the differences are between -0.003～0.003 while a search range of the axial strain is about 0.9～1.1. This shows that the GPU implementation achieves nearly the same result as the CPU implementation. The minor differences may be caused by numerical computing differences accumulated during further computation.

[0137] As for memory usage, the GPU implementation uses less than 4 GB of device memory and less than 12GB of main memory. As for compatibility, the experiment uses the Nvidia Fourier Transform library cuFFT. All other included libraries are common C++ libraries or helper function provided by Nvidia. The program can be ported to other platforms or systems with corresponding compilers.

[0138] Fig. 21 depicts an exemplary computing device 2100, hereinafter interchangeably referred to as a computer system 2100. The following description of the computing device 2100 is provided by way of example only and is not intended to be limiting.

[0139] As shown in Fig. 21, the example computing device 2100 includes a processor 2102 for executing software routines. Although a single processor is shown for the sake of clarity, the computing device 2100 may also include a multi-processor system. The processor 2102 may include the CPU and / or the GPU. The processor 2102 is connected to a communication infrastructure 2104 for communication with other components of the computing device 2100. The communication infrastructure 2104 may include, for example, a communications bus, cross-bar, or network.

[0140] The computing device 2100 further includes a main memory 2106, such as a random access memory (RAM) , and a secondary memory 2108. The secondary memory 2108 may include, for example, a hard disk drive 2110 and / or a removable storage drive 2112, which may include a floppy disk drive, a magnetic tape drive, an optical disk drive, or the like. The removable storage drive 2112 reads from and / or writes to a removable storage unit 2114 in a well-known manner. The removable storage unit 2114 may include a floppy disk, magnetic tape, optical disk, or the like, which is read by and written to by removable storage drive 2112. As will be appreciated by persons skilled in the relevant art (s) , the removable storage unit 2114 includes a computer readable storage medium having stored therein computer executable program code instructions and / or data.

[0141] In an alternative implementation, the secondary memory 2108 may additionally or alternatively include other similar means for allowing computer programs or other instructions to be loaded into the computing device 2100. Such means can include, for example, a removable storage unit 2116 and an interface 2118. Examples of a removable storage unit 2116 and interface 2118 include a program cartridge and cartridge interface (such as that found in video game console devices) , a removable memory chip (such as an EPROM or PROM) and associated socket, and other removable storage units 2116 and interfaces 2118 which allow software and data to be transferred from the removable storage unit 2116 to the computer system 2100.

[0142] The computing device 2100 also includes at least one communication interface 2120. The communication interface 2120 allows software and data to be transferred between computing device 2100 and external devices via a communication path 2122. In various embodiments of the inventions, the communication interface 2120 permits data to be transferred between the computing device 2100 and a data communication network, such as a public data or private data communication network. The communication interface 2120 may be used to exchange data between different computing devices 2100 which such computing devices 2100 form part an interconnected computer network. Examples of a communication interface 2120 can include a modem, a network interface (such as an Ethernet card) , a communication port, an antenna with associated circuitry and the like. The communication interface 2120 may be wired or may be wireless. Software and data transferred via the communication interface 2120 are in the form of signals which can be electronic, electromagnetic, optical or other signals capable of being received by communication interface 2120. These signals are provided to the communication interface via the communication path 2122.

[0143] As shown in Fig. 21, the computing device 2100 further includes a display interface 2124 which performs operations for rendering images to an associated display 2126 and an audio interface 2128 for performing operations for playing audio content via associated speaker (s) 2130.

[0144] As used herein, the term "computer program product" may refer, in part, to removable storage unit 2114, removable storage unit 2116, a hard disk installed in hard disk drive 2110, or a carrier wave carrying software over communication path 2122 (wireless link or cable) to communication interface 2120. Computer readable storage media refers to any non-transitory tangible storage medium that provides recorded instructions and / or data to the computing device 2100 for execution and / or processing. Examples of such storage media include floppy disks, magnetic tape, CD-ROM, DVD, Blu-rayTM Disc, a hard disk drive, a ROM or integrated circuit, USB memory, a magneto-optical disk, or a computer readable card such as a PCMCIA card and the like, whether or not such devices are internal or external of the computing device 2100. Examples of transitory or non-tangible computer readable transmission media that may also participate in the provision of software, application programs, instructions and / or data to the computing device 2100 include radio or infra-red transmission channels as well as a network connection to another computer or networked device, and the Internet or Intranets including e-mail transmissions and information recorded on Websites and the like.

[0145] The computer programs (also called computer program code) are stored in main memory 2106 and / or secondary memory 2108. Computer programs can also be received via the communication interface 2120. Such computer programs, when executed, enable the computing device 2100 to perform one or more features of embodiments discussed herein. In various embodiments, the computer programs, when executed, enable the processor 2102 to perform features of the above-described embodiments. Accordingly, such computer programs represent controllers of the computer system 2100.

[0146] Software may be stored in a computer program product and loaded into the computing device 2100 using the removable storage drive 2112, the hard disk drive 2110, or the interface 2118. Alternatively, the computer program product may be downloaded to the computer system 2100 over the communications path 2122. The software, when executed by the processor 2102, causes the computing device 2100 to perform functions of embodiments described herein.

[0147] It is to be understood that the embodiment of Fig. 21 is presented merely by way of example. Therefore, in some embodiments one or more features of the computing device 2100 may be omitted. Also, in some embodiments, one or more features of the computing device 2100 may be combined together. Additionally, in some embodiments, one or more features of the computing device 2100 may be split into one or more component parts.

[0148] It will be appreciated that the elements illustrated in Fig. 21 function to provide means for performing the various functions and operations of the servers as described in the above embodiments.

[0149] In an implementation, a server may be generally described as a physical device comprising at least one processor and at least one memory including computer program code. The at least one memory and the computer program code are configured to, with the at least one processor, cause the physical device to perform the requisite operations.

[0150] When the computing device 2100 is configured to realise the system 700 for performing an accelerated coupled filtering method for tissue deformation analysis, the system 700 can have a non-transitory computer readable medium having stored thereon an application which when executed causes the system 700 to perform the accelerated coupled filtering method 800 described above and / or any other methods described herein.

[0151] It will be appreciated by a person skilled in the art that numerous variations and / or modifications may be made to the present invention as shown in the specific embodiments without departing from the spirit or scope of the invention as broadly described. The present embodiments are, therefore, to be considered in all respects to be illustrative and not restrictive.

Claims

1.An accelerated coupled filtering method for tissue deformation analysis, comprising:initializing, by a first processing device, a plurality of first motion parameters;applying, by a second processing device, a first filter on a pre-deformation image of a tissue to obtain a filtered pre-deformation image, wherein the second processing device is different from the first processing device;applying, by the second processing device, a second filter on a post-deformation image of the tissue to obtain a filtered post-deformation image, wherein the filtered pre-deformation image and the filtered post-deformation image are correlated by a first motion matrix comprising the plurality of first motion parameters, and wherein the plurality of first motion parameters comprise at least three fundamental first motion parameters, each of the three fundamental first motion parameters represents movement of the tissue along an axial direction relative to the axial direction, an elevational direction and a lateral direction, respectively;estimating, by the second processing device, a respective value for each of the plurality of first motion parameters, wherein each of the estimated respective values represents a difference between the pre-deformation and post-deformation filtered images; andin response to determining that at least one of the estimated respective values meets at least one predefined criterion, updating the at least one of the estimated respective values as an optimal value for at least one corresponding first motion parameter of the plurality of first motion parameters.2.The method of claim 1, wherein the first processing device is a Central Processing Unit (CPU) and the second processing device is a Graphical Processing Unit (GPU) .3.The method of claim 2, further comprising executing, by the GPU, at least one kernel function associated with multi-threading on a plurality of computing cores.4.The method of claim 1, further comprising:detecting, by the second processing device, one or more envelopes present in the filtered pre-deformation image and filtered post-deformation image to obtain a filtered B-mode pre-deformation image and a filtered B-mode post-deformation image.5.The method of claim 1, wherein applying the first filter and the second filter comprises convolving a first point spread function and second point spread function of an ultrasound system with the pre-deformation image and the post-deformation image, respectively, and wherein the first point spread function is a modified version of the second point spread function.6.The method of claim 1, wherein, after applying the second filter, the method further comprises, spatially transforming, by the second processing device, the post-deformation image based at least on a second motion matrix, and wherein the first motion matrix is a modified version of the second motion matrix.7.The method of claim 6, wherein the filtered pre-deformation image and filtered post-deformation image are 3-dimensional (3D) images, and wherein estimating the respective values further comprises:(A) flattening, by the second processing device, a first and second plurality of voxels surrounding each of one or more pre-deformation points of interest in the filtered pre-deformation image and each of one or more post-deformation points of interest in the filtered post-deformation image, respectively, wherein the flattening is performed along the elevational direction of each of the first and second plurality of voxels, and wherein the flattening reduces the plurality of first and second motion parameters.8.The method of claim 7, wherein estimating the respective values is performed in a plurality of iterations, and wherein:in a first iteration of the plurality of iterations:the plurality of first motion parameters comprise two of the three fundamental first motion parameters representing movement of the tissue along the axial direction relative to the axial direction and the lateral direction, respectively, and wherein the remaining fundamental first motion parameter is removed by the flattening, andbetween a second iteration to a final iteration of the plurality of iterations:the plurality of first motion parameters comprise the two of the three fundamental first motion parameters, two lateral first motion parameters representing movement of the tissue along the lateral direction relative to the lateral direction and the axial direction, respectively, and two elevational first motion parameters representing movement of the tissue along the elevational direction relative to the lateral direction and the axial direction, respectively.9.The method of claim 8, further comprising executing, by the second processing device, one or more Fast Fourier Transform batch plans, wherein the one or more Fast Fourier Transform batch plans are configured to:perform Fourier Transform on a plurality of images in the first iteration; orperform Fourier Transform on the first and second plurality of voxels surrounding each of the one or more pre-deformation points of interest in the filtered pre-deformation image and each of the one or more post-deformation points of interest in the filtered post-deformation image, concurrently.10.The method of claim 8, wherein the second motion matrix comprises a plurality of second motion parameters, and wherein estimating the respective values further comprises:(B) searching, by the second processing device, a value of each of the plurality of first motion parameters for each of the one or more pre-deformation points of interest in the filtered pre-deformation image;(C) searching, by the second processing device, a value of each of the plurality of second motion parameters for each of the one or more post-deformation points of interest in the filtered post-deformation image;(D) generating, by the second processing device, one or more pre-deformation voxels based on the searched value of each of the plurality of first motion parameters;(E) generating, by the second processing device, one or more post-deformation voxels based on the searched value of each of the plurality of second motion parameters;(F) calculating, by the second processing device, a similarity metric based on one or more first voxels surrounding the one or more pre-deformation voxels and one or more second voxels surrounding the one or more post-deformation voxels, wherein the similarity metric defines a similarity between a target pre-deformation point of interest and a post-deformation point of interest corresponding to the target pre-deformation point of interest; and(G) determining, by the second processing device, if the target pre-deformation point of interest matches the corresponding post-deformation point of interest based on the calculated similarity metric.11.The method of claim 10, wherein the plurality of second motion parameters comprise at least two fundamental second motion parameters, two lateral second motion parameters and two elevational motion parameters, and wherein estimating the respective values comprises:in the first iteration:determining, by the second processing device, a value of each of the two of the three fundamental first motion parameters based on steps (A) to (G) ;between the second iteration to a third final iteration:determining, by the second processing device, an updated value of each of the two of the three fundamental first motion parameters based on steps (A) to (G) ; andassigning, by the second processing device, a value of each of the two lateral second motion parameters and the two elevational second motion parameters searched in step (C) in a previous iteration as a respective value of each of the two lateral first motion parameters and the two elevational first motion parameters; and in a final two iterations:assigning, by the second processing device, the determined updated value of each of the two of the three fundamental first motion parameters in the third final iteration as an optimal value of each of the two of the three fundamental first motion parameters; andassigning, by the second processing device, an updated value of each of the two lateral second motion parameters and the two elevational second motion parameters searched in step (C) in the third final iteration as an optimal value of each of the two lateral first motion parameters and the two elevational first motion parameters.12.The method of claim 11, wherein the similarity metric comprises any one of a Normalized Correlation Coefficient, Sum of Absolute Differences and Sum of Squared Differences.13.The method of claim 1, further comprising:imposing, by the second processing device, an upper limit and a lower limit on a gradient of a displacement range of the fundamental first motion parameters along the axial direction, wherein values falling outside the displacement range are determined to be a wrong value of the fundamental first motion parameters;in response to determining that one or more voxels comprise the wrong value of the fundamental first motion parameters:determining, by the second processing device, a corrected value of the fundamental first motion parameters of the one or more voxels, wherein the corrected value is determined by interpolating the wrong value based on a predetermined correct value of the fundamental first motion parameters of voxels surrounding the one or more voxels.14.A system for performing an accelerated coupled filtering method for tissue deformation analysis, the system comprising a first processing device and a second processing device, wherein the second processing device is different from the first processing device, andwherein the first processing device is configured to:initialize a plurality of first motion parameters; andwherein the second processing device is configured to:apply a first filter on a pre-deformation image of a tissue to obtain a filtered pre-deformation image;apply a second filter on a post-deformation image of the tissue to obtain a filtered post-deformation image, wherein the filtered pre-deformation image and the filtered post-deformation image are correlated by a first motion matrix comprising the plurality of first motion parameters, and wherein the plurality of first motion parameters comprise at least three fundamental first motion parameters, each of the three fundamental first motion parameters represents movement of the tissue along an axial direction relative to the axial direction, an elevational direction and a lateral direction, respectively;estimate a respective value for each of the plurality of first motion parameters, wherein each of the estimated respective values represent a difference between the pre-deformation and post-deformation filtered images; andin response to determining that at least one of the estimated respective values meets at least one predefined criterion, update the at least one of the estimated respective values as an optimal value for at least one corresponding first motion parameter of the plurality of first motion parameters.15.The system of claim 14, wherein the first processing device is a Central Processing Unit (CPU) and the second processing device is a Graphical Processing Unit (GPU) .16.The system of claim 15, wherein the GPU is configured to execute at least one kernel function associated with multi-threading on a plurality of computing cores.17.The system of claim 14, wherein the second processing device is further configured to:detect one or more envelopes present in the filtered pre-deformation image and filtered post-deformation image to obtain a filtered B-mode pre-deformation image and a filtered B-mode post-deformation image.18.The system of claim 14, wherein to apply the first filter and the second filter, the second processing device is configured to convolve a first point spread function and second point spread function of an ultrasound system with the pre-deformation image and the post-deformation image, respectively, and wherein the first point spread function is a modified version of the second point spread function.19.The system of claim 14, wherein, after applying the second filter, the second processing device is further configured to spatially transform the post-deformation image based at least on a second motion matrix, and wherein the first motion matrix is a modified version of the second motion matrix.20.The system of claim 19, wherein the filtered pre-deformation image and filtered post-deformation image are 3-dimensional (3D) images, and wherein to estimate the respective values, the second processing device is configured to:(A) flatten a first and second plurality of voxels surrounding each of one or more pre-deformation points of interest in the filtered pre-deformation image and each of one or more post-deformation points of interest in the filtered post-deformation image, respectively, wherein the flattening is performed along the elevational direction of each of the first and second plurality of voxels, and wherein the flattening reduces the plurality of first and second motion parameters.21.The system of claim 20, wherein estimating the respective values is performed in a plurality of iterations, and wherein:in a first iteration of the plurality of iterations:the plurality of first motion parameters comprise two of the three fundamental first motion parameters representing movement of the tissue along the axial direction relative to the axial direction and the lateral direction, respectively, and wherein the remaining fundamental first motion parameter is removed by the flattening, andbetween a second iteration to a final iteration of the plurality of iterations:the plurality of first motion parameters comprise the two of the three fundamental first motion parameters, two lateral first motion parameters representing movement of the tissue along the lateral direction relative to the lateral direction and the axial direction, respectively, and two elevational first motion parameters representing movement of the tissue along the elevational direction relative to the lateral direction and the axial direction, respectively.22.The system of claim 21, wherein the second processing device is further configured to execute one or more Fast Fourier Transform batch plans, and wherein the one or more Fast Fourier Transform batch plans are configured to:perform Fourier Transform on a plurality of images in the first iteration; or perform Fourier Transform on the first and second plurality of voxels surrounding each of the one or more pre-deformation points of interest in the filtered pre-deformation image and each of the one or more post-deformation points of interest in the filtered post-deformation image, concurrently.23.The system of claim 21, wherein the second motion matrix comprises a plurality of second motion parameters, and wherein to estimate the respective values, the second processing device is further configured to:(B) search a value of each of the plurality of first motion parameters for each of the one or more pre-deformation points of interest in the filtered pre-deformation image;(C) search a value of each of the plurality of second motion parameters for each of the one or more post-deformation points of interest in the filtered post-deformation image;(D) generate one or more pre-deformation voxels based on the searched value of each of the plurality of first motion parameters;(E) generate one or more post-deformation voxels based on the searched value of each of the plurality of second motion parameters;(F) calculate a similarity metric based on one or more first voxels surrounding the one or more pre-deformation voxels and one or more second voxels surrounding the one or more post-deformation voxels, wherein the similarity metric defines a similarity between a target pre-deformation point of interest and a post-deformation point of interest corresponding to the target pre-deformation point of interest; and(G) determine if the target pre-deformation point of interest matches the corresponding post-deformation point of interest based on the calculated similarity metric.24.The system of claim 23, wherein the plurality of second motion parameters comprise at least two fundamental second motion parameters, two lateral second motion parameters and two elevational motion parameters, and wherein to estimate the respective values, the processing device is further configured to:in the first iteration:determine a value of each of the two of the three fundamental first motion parameters based on steps (A) to (G) ;between the second iteration to a third final iteration:determine an updated value of each of the two of the three fundamental first motion parameters based on steps (A) to (G) ; andassign a value of each of the two lateral second motion parameters and the two elevational second motion parameters searched in step (C) in a previous iteration as a respective value of each of the two lateral first motion parameters and the two elevational first motion parameters; andin a final two iterations:assign the determined updated value of each of the two of the three fundamental first motion parameters in the third final iteration as an optimal value of each of the two of the three fundamental first motion parameters; andassign an updated value of each of the two lateral second motion parameters and the two elevational second motion parameters searched in step (C) in the third final iteration as an optimal value of each of the two lateral first motion parameters and the two elevational first motion parameters.25.The system of claim 24, wherein the similarity metric comprises any one of a Normalized Correlation Coefficient, Sum of Absolute Differences and Sum of Squared Differences.26.The system of claim 14, wherein the second processing device is further configured to:impose an upper limit and a lower limit on a gradient of a displacement range of the fundamental first motion parameters along the axial direction, wherein values falling outside the displacement range are determined to be a wrong value of the fundamental first motion parameters;in response to determining that one or more voxels comprise the wrong value of the fundamental first motion parameters:determine a corrected value of the fundamental first motion parameters of the one or more voxels, wherein the corrected value is determined by interpolating the wrong value based on a predetermined correct value of the fundamental first motion parameters of voxels surrounding the one or more voxels.

Images