Accelerated coupling filtering method and system for tissue deformation analysis

Abstract

There is provided an accelerated coupled filtering method for tissue deformation analysis, comprising the steps of applying, by a processing device, a first filter to a pre-deformation image of tissue to obtain a filtered pre-deformation image, and applying, by the processing device, a second filter to a post-deformation image of the tissue to obtain a filtered post-deformation image. The filtered pre-deformation image and the filtered post-deformation image may be associated by a first motion matrix including a plurality of first motion parameters. The plurality of first motion parameters may include at least three substantially first motion parameters, each of which may represent movement of the tissue along an axial direction relative to an axial direction, a vertical direction, and a lateral direction, respectively. The method further comprises the step of estimating, by the processing device, a respective value for each of the plurality of first motion parameters. Each of the estimated respective values may represent a difference between the pre-deformation filtered image and the post-deformation filtered image.

Description

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims priority to U.S. Provisional Patent Application No. 63 / 567,427, filed on March 20, 2024, which is incorporated herein by reference in its entirety. Technical Field

[0003] The present invention relates broadly, but not exclusively, to acceleration coupled filtering methods and systems for tissue deformation analysis. Background Art

[0004] Ultrasound image-based tissue deformation analysis can be used to measure tissue stiffness and reveal useful information for clinical diagnosis. This analysis is typically performed by comparing images taken before and after tissue deformation. However, feature motion decorrelation caused by the ultrasound imaging procedure can greatly hinder accuracy. To address this issue, a coupled filtering method was proposed and implemented to analytically compensate for feature motion decorrelation. Although the coupled filtering method can achieve much higher accuracy than other existing methods, this method is computationally expensive. Therefore, the implementation of this 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 would require billions of times more operations than the 2D implementation, making the 3D implementation impractical for real-world applications.

[0005] Therefore, there is a need to provide an acceleration coupled filtering method and system for tissue deformation analysis. Summary of the Invention

[0006] According to a first aspect of the present invention, an accelerated coupled filtering method for tissue deformation analysis is provided. The method comprises: applying a first filter to a pre-deformation image of the tissue by a processing device to obtain a filtered pre-deformation image; applying a second filter to a post-deformation image of the tissue by the processing device 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 a plurality of first motion parameters, wherein the plurality of first motion parameters comprises at least three basic first motion parameters, each of the three basic first motion parameters representing movement of the tissue along an axial direction relative to an axial direction, a vertical direction, and a transverse direction, respectively; estimating, by the processing device, a corresponding value for each of the plurality of first motion parameters, each of the estimated corresponding values ​​representing a difference between the pre-deformation filtered image and the post-deformation filtered image; and, in response to determining that at least one of the estimated corresponding values ​​satisfies at least one predefined criterion, updating at least one of the estimated corresponding values ​​to an optimal value of at least one corresponding first motion parameter of the plurality of first motion parameters.

[0007] In an embodiment of the present disclosure, applying the first filter and the second filter may include convolving first and second point spread functions of the ultrasound system with the pre-deformation image and the deformed image, respectively, and the first point spread function may be a modified version of the second point spread function.

[0008] In an embodiment of the present disclosure, the method may include, after applying the second filter, a step of spatially transforming the warped image by the processing device based at least on the second motion matrix. The first motion matrix may be a modified version of the second motion matrix.

[0009] In an embodiment of the present disclosure, the method may include a step of detecting, by a processing device, one or more envelopes present in a filtered pre-deformation image and a filtered post-deformation image to obtain a filtered B-mode pre-deformation image and a filtered B-mode post-deformation image.

[0010] In an embodiment of the present disclosure, the filtered pre-deformation image and the filtered post-deformation image may be three-dimensional (3D) images, and estimating the corresponding value may include: (A) flattening, by a processing device, a first plurality of voxels and a second plurality of voxels around each of one or more pre-deformation interest points in the filtered pre-deformation image and each of one or more post-deformation interest points in the filtered post-deformation image, respectively. The flattening may be performed along a vertical direction of each of the first plurality of voxels and the second plurality of voxels, and the flattening may reduce a plurality of first motion parameters and a plurality of second motion parameters.

[0011] In an embodiment of the present disclosure, estimating the corresponding value can be performed in multiple iterations. In the first iteration of the multiple iterations, the multiple first motion parameters may include two basic first motion parameters of the three basic first motion parameters that respectively represent the movement of the tissue along the axial direction relative to the axial direction and the transverse direction. The remaining basic first motion parameters can be removed by flattening. Between the second iteration and the final iteration of the multiple iterations, the multiple first motion parameters may include two basic first motion parameters of the three basic first motion parameters, two transverse first motion parameters that respectively represent the movement of the tissue along the transverse direction relative to the transverse direction and the axial direction, and two vertical first motion parameters that respectively represent the movement of the tissue along the vertical direction relative to the transverse direction and the axial direction.

[0012] In an embodiment of the present disclosure, the second motion matrix may include a plurality of second motion parameters. Estimating the corresponding value may also include: (B) searching, by the processing device, for a value of each of the plurality of first motion parameters for each of the one or more pre-deformed interest points in the filtered pre-deformed image; (C) searching, by the processing device, for a value of each of the plurality of second motion parameters for each of the one or more post-deformed interest points in the filtered post-deformed image; (D) generating, by the processing device, one or more pre-deformed voxels based on the search value of each of the plurality of first motion parameters; (E) generating, by the processing device, one or more post-deformed voxels based on the search value of each of the plurality of second motion parameters; (F) calculating, by the processing device, a similarity index based on one or more first voxels surrounding the one or more pre-deformed voxels and one or more second voxels surrounding the one or more post-deformed voxels, the similarity index defining a similarity between the target pre-deformed interest point and the post-deformed interest point corresponding to the target pre-deformed interest point; and (G) determining, by the processing device, whether the target pre-deformed interest point matches the corresponding post-deformed interest point based on the calculated similarity index.

[0013] In an embodiment of the present disclosure, the plurality of second motion parameters may include at least two basic second motion parameters, two lateral second motion parameters, and two vertical motion parameters. Estimating the corresponding values ​​may also include: in the first iteration, determining, by the processing device, the value of each of two basic first motion parameters among the three basic first motion parameters based on steps (A) to (G); between the second iteration and the third final iteration, determining, by the processing device, an updated value of each of two basic first motion parameters among the three basic first motion parameters based on steps (A) to (G), and assigning, by the processing device, the value of each of the two horizontal second motion parameters and the two vertical second motion parameters searched in step (C) in the previous iteration as the corresponding value of each of the two horizontal first motion parameters and the two vertical first motion parameters; and in the final two iterations, assigning, by the processing device, the determined updated value of each of the two basic first motion parameters among the three basic first motion parameters in the third final iteration as the optimal value of each of the two basic first motion parameters among the three basic first motion parameters, and assigning, by the processing device, the updated value of each of the two horizontal second motion parameters and the two vertical second motion parameters searched in step (C) in the third final iteration as the optimal value of each of the two horizontal first motion parameters and the two vertical first motion parameters.

[0014] In an embodiment of the present disclosure, the similarity index may include any one of a normalized correlation coefficient, a sum of absolute differences, and a sum of squared differences.

[0015] In an embodiment of the present disclosure, the method may further include: a processing device imposing an upper limit and a lower limit on the gradient of the displacement range of the basic first motion parameter along the axial direction, wherein a value falling outside the displacement range may be determined as an erroneous value of the basic first motion parameter; in response to determining that one or more voxels include an erroneous value of the basic first motion parameter, the processing device determines a correction value of the basic first motion parameter of the one or more voxels, and the correction value may be determined by interpolating the erroneous value based on a predetermined correct value of the basic first motion parameter of the voxels surrounding the one or more voxels.

[0016] According to a second aspect of the present invention, a system for an accelerated coupled filtering method for performing deformation analysis of tissue is provided, the system including a processing device configured to: apply a first filter to a pre-deformation image of the tissue to obtain a filtered pre-deformation image; apply a second filter to the post-deformation image of the tissue to obtain a filtered post-deformation image, the filtered pre-deformation image and the filtered post-deformation image being related by a first motion matrix including a plurality of first motion parameters, and wherein the plurality of first motion parameters include at least three basic first motion parameters, each of the three basic first motion parameters representing the movement of the tissue along the axial direction relative to the axial direction, the vertical direction and the transverse direction, respectively; estimate a corresponding value of each of the plurality of first motion parameters, each of the estimated corresponding values ​​representing a difference between the pre-deformation filtered image and the post-deformation filtered image; and in response to determining that at least one of the estimated corresponding values ​​meets at least one predefined criterion, update at least one of the estimated corresponding values ​​to an optimal value of at least one first motion parameter corresponding to the plurality of first motion parameters.

[0017] In an embodiment of the present disclosure, applying the first filter and the second filter may include convolving a first point spread function and a second point spread function of the ultrasound system with the pre-deformation image and the deformed image, respectively, and the first point spread function may be a modified version of the second point spread function.

[0018] In an embodiment of the present disclosure, the system may be configured to: after applying the second filter, spatially transform the warped image based on at least the second motion matrix.The first motion matrix may be a modified version of the second motion matrix.

[0019] In an embodiment of the present disclosure, the processing device may be configured to detect one or more envelopes present in the filtered pre-warping image and the filtered post-warping image to obtain the filtered B-mode pre-warping image and the filtered B-mode post-warping image.

[0020] In an embodiment of the present disclosure, the filtered pre-deformation image and the filtered post-deformation image may be three-dimensional (3D) images, and to estimate the corresponding value, the processing device may be configured to: (A) flatten a first plurality of voxels and a second plurality of voxels around 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 may be performed along a vertical direction of each of the first plurality of voxels and the second plurality of voxels, and the flattening may reduce a plurality of first motion parameters and a plurality of second motion parameters.

[0021] In an embodiment of the present disclosure, estimating the corresponding value can be performed in multiple iterations. In the first iteration of the multiple iterations, the multiple first motion parameters may include two basic first motion parameters of the three basic first motion parameters that respectively represent the movement of the tissue along the axial direction relative to the axial direction and the transverse direction. The remaining basic first motion parameters can be removed by flattening. Between the second iteration and the final iteration of the multiple iterations, the multiple first motion parameters may include two basic first motion parameters of the three basic first motion parameters, two transverse first motion parameters that respectively represent the movement of the tissue along the transverse direction relative to the transverse direction and the axial direction, and two vertical first motion parameters that respectively represent the movement of the tissue along the vertical direction relative to the transverse direction and the axial direction.

[0022] In an embodiment of the present disclosure, the second motion matrix may include a plurality of second motion parameters. In addition, to estimate the corresponding value, the processing device may be further configured to: (B) search for a value of each of the plurality of first motion parameters for each of the one or more pre-deformed interest points in the filtered pre-deformed image; (C) search for a value of each of the plurality of second motion parameters for each of the one or more post-deformed interest points in the filtered post-deformed image; (D) generate one or more pre-deformed voxels based on the search value of each of the plurality of first motion parameters; (E) generate one or more post-deformed voxels based on the search value of each of the plurality of second motion parameters; (F) calculate a similarity index based on one or more first voxels surrounding the one or more pre-deformed voxels and one or more second voxels surrounding the one or more post-deformed voxels, the similarity index being able to define a similarity between the target pre-deformed interest point and the post-deformed interest point corresponding to the target pre-deformed interest point; and (G) determine whether the target pre-deformed interest point matches the corresponding post-deformed interest point based on the calculated similarity index.

[0023] In an embodiment of the present disclosure, the plurality of second motion parameters may include at least two basic second motion parameters, two lateral second motion parameters, and two vertical motion parameters. In order to estimate the implementation corresponding values, the processing device can also be configured to: in the first iteration, determine the value of each of two basic first motion parameters among the three basic first motion parameters based on steps (A) to (G); between the second iteration and the third final iteration, determine the updated value of each of two basic first motion parameters among the three basic first motion parameters based on steps (A) to (G), and assign the value of each of the two horizontal second motion parameters and the two vertical second motion parameters searched in step (C) in the previous iteration as the corresponding value of each of the two horizontal first motion parameters and the two vertical first motion parameters; and in the final two iterations, assign the determined updated value of each of the two basic first motion parameters among the three basic first motion parameters in the third final iteration as the optimal value of each of the two basic first motion parameters among the three basic first motion parameters, and assign the updated value of each of the two horizontal second motion parameters and the two vertical second motion parameters searched in step (C) in the third final iteration as the optimal value of each of the two horizontal first motion parameters and the two vertical first motion parameters.

[0024] In an embodiment of the present disclosure, the similarity index may include any one of a normalized correlation coefficient, a sum of absolute differences, and a sum of squared differences.

[0025] In an embodiment of the present disclosure, the processing device can be configured to: impose an upper limit and a lower limit on the gradient of the displacement range of the basic first motion parameter along the axial direction, and the value falling outside the displacement range can be determined as an erroneous value of the basic first motion parameter; in response to determining that one or more voxels include an erroneous value of the basic first motion parameter, determine a correction value of the basic first motion parameter of the one or more voxels, and the correction value can be determined by interpolating the erroneous value based on a predetermined correct value of the basic first motion parameter of the voxels surrounding the one or more voxels. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Those skilled in the art will better understand and readily appreciate the various embodiments of the present invention through the following written description, which is given by way of example only, taken in conjunction with the accompanying drawings, in which:

[0027] Figure 1 A flow chart illustrating a conventional tissue deformation analysis process is shown.

[0028] Figure 2A The left side of shows an intuitive diagram for illustrating the coarse-to-fine search strategy, and Figure 2A The right side of FIG. 1 shows a flow chart illustrating a coarse-to-fine search strategy for a conventional “baseline” coupling filtering method.

[0029] Figure 2B A flow chart illustrating the detailed search process of a conventional "baseline" coupling filtering method is shown.

[0030] Figure 3 A table listing parameters used in FastCF according to an embodiment is shown.

[0031] Figure 4 Shown are flowcharts for illustrating (a) a conventional process of converting a radio frequency (RF) image into a B-mode image and (b) a process of converting an RF image into a B-mode image according to an embodiment.

[0032] Figure 5 An example scenario illustrating a post-processing step according to an embodiment.

[0033] Figure 6 Shown is a flowchart illustrating an overview of an accelerated coupled filtering method ("FastCF") according to an embodiment.

[0034] Figure 7 A schematic diagram illustrating an exemplary system for performing an accelerated coupled filtering method for tissue deformation analysis according to an embodiment.

[0035] Figure 8 Shown is a flow chart illustrating an exemplary workflow of an accelerated coupled filtering method for tissue deformation analysis, according to an embodiment.

[0036] Figure 9 An exemplary model of elastic tissue according to an embodiment is shown.

[0037] Figure 10 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 is shown, according to an embodiment.

[0038] Figure 11 A table listing parameters for time complexity and space complexity analysis of 2D and 3D implementations of various coupled filtering methods is shown, according to an embodiment.

[0039] Figure 12 Another table is shown listing detailed time complexity analysis of 3D implementations of various coupled filtering methods, according to an embodiment.

[0040] Figure 13 A table listing the square root of mean square error (SRMSE) and computation time of FastCF performed on a 3D simulation dataset is shown, according to an embodiment.

[0041] Figure 14Shown is a color map (in grayscale) illustrating results obtained from executing FastCF on a 3D simulated dataset, according to an embodiment.

[0042] Figure 15 A table listing the computation times of different functions used in the FastCF process on a 3D simulation dataset is shown, according to an embodiment.

[0043] Figure 16 Shown is a color map (in grayscale) illustrating the results obtained from performing FastCF on a 3D real dataset, according to an embodiment.

[0044] Figure 17 Color maps (in grayscale) showing (a) a baseline method and (b) FastCF performed on a 2D simulated dataset, according to an embodiment.

[0045] Figure 18 Shown is a graph illustrating the Normalized Correlation Coefficient (NCC) of images based on different motion parameters used in FastCF, according to an embodiment.

[0046] Figure 19 A table listing the SRMSE of FastCF under different setting options according to an embodiment is shown.

[0047] Figure 20 A schematic diagram illustrating an example of a computing device for implementing a system for performing an accelerated coupled filtering method for tissue deformation analysis, according to an embodiment. DETAILED DESCRIPTION

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

[0049] Some portions of the description that follows are presented, either explicitly or implicitly, 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 most effectively convey the substance of their work to others skilled in the art. An algorithm is here, and generally, considered to be a self-consistent sequence of steps leading to a desired result. Such steps are those requiring physical manipulation of physical quantities, such as electrical, magnetic, or optical signals capable of being stored, transferred, combined, compared, and otherwise manipulated.

[0050] Unless otherwise specifically noted, and as will become apparent from the following text, it will be understood that throughout this specification, discussions utilizing terms such as "scan," "calculate," "determine," "replace," "generate," "initialize," "output," and the like refer to actions and processes of a computer system or similar electronic device that manipulate and transform 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.

[0051] This specification also discloses a device for performing the operations of the method. Such a device may be specially constructed for the desired purpose, or may include 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 device. Various machines can be used with the programs taught herein. Alternatively, it may be appropriate to construct a more specialized device to perform the required method steps. The structure of a conventional computer will become apparent from the following description.

[0052] In addition, this specification also implicitly discloses a computer program, because it is obvious to those skilled in the art that each step of the method described herein can be implemented by computer code. The computer program is not intended to be limited to any specific programming language and its implementation method. It is understood that the teachings of the disclosure contained herein can be implemented using a variety of programming languages ​​and their encoding. In addition, the implementation of the computer program is not intended to be limited to any specific control flow. Without departing from the spirit or scope of the present invention, there are many other variations of the computer program, which use different control flows.

[0053] In addition, one or more steps of the computer program can be performed in parallel rather than sequentially. Such a computer program can be stored on any computer-readable medium. Computer-readable media can include storage devices, such as magnetic disks or optical disks, memory chips or other storage devices suitable for being connected to a computer interface. Computer-readable media can also include hard-wired media (such as exemplified in an Internet system) or wireless media (such as exemplified in a GSM, GPRS, 3G or 4G mobile phone system) and other wireless systems (such as Bluetooth, ZigBee, Wi-Fi). The computer program, when loaded on such a computer and executed, effectively generates a device for implementing the steps of the preferred method.

[0054] The present invention may also be implemented as a hardware module. More specifically, in the hardware sense, a module is a functional hardware unit designed to be used in conjunction with other components or modules. For example, a module may be implemented using discrete electronic components, or a module may form part of an entire electronic circuit, such as an application-specific integrated circuit (ASIC) or a field-programmable gate array (FPGA). Many other possibilities exist. Those skilled in the art will appreciate that the system may also be implemented as a combination of hardware and software modules.

[0055] In the following description, the term "module" may refer to a software element, a hardware element, or a combination of both.

[0056] An application programming interface (API) enables software and applications to communicate with each other. It is a software-to-software interface that allows independent parties to communicate with each other without any prior user knowledge or intervention. Generally speaking, it is a set of well-defined communication methods between various software components.

[0057] This application uses the term "configured to" when referring to systems, devices, and computer program components. For a system of one or more computers configured to perform a particular operation or action, it means that the system has installed on it software, firmware, hardware, or a combination thereof that causes the system to perform those operations or actions in operation. For one or more computer programs configured to perform a particular operation or action, it means that the implementation of one or more programs includes instructions that, when executed by a data processing device, cause the device to perform those operations or actions. For a dedicated logic circuit configured to perform a particular operation or action, it means that the circuit has electronic logic that performs those operations or actions.

[0058] As used herein, the term "processing device" refers to any hardware or system configured to perform a computing task.

[0059] For example, malignant tumors, cirrhosis of the liver, and dead heart tissue are significantly stiffer than healthy tissue. Therefore, stiffness can reveal pathology, leading doctors to perform tissue deformation analysis to track deformation caused by internal or external forces in order to infer stiffness distribution for diagnosis. Figure 1A flowchart illustrating a typical process of tissue deformation analysis is shown. For accurate tissue deformation analysis, a coupled filtering method is proposed, which applies two different but coupled filters to images taken before and after tissue deformation to compensate for both complex motion (i.e., feature motion decorrelation) and echo interference. The coupled filtering method can achieve much higher accuracy than other correlation methods, especially for large tissue deformations. However, this method may require re-filtering images with different motions by exhaustively searching all possible discrete motions. Therefore, such a heavy computational load makes it feasible to implement it only on 2D images in practical use. Analyzing a pair of images of size 201×1001 in MATLAB on a CPU may take about 45 hours, while it takes about 20 minutes on an FPGA.

[0060] In addition, since the tissue moves 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 of all possible affine motions, the 3D implementation of the coupled filtering method is estimated to perform billions of more calculations than the 2D implementation. Therefore, the present disclosure seeks to improve the algorithmic complexity of the coupled filtering method. In particular, the present disclosure provides a method (hereinafter referred to as "FastCF") that can accelerate the coupled filtering method while maintaining similar accuracy. The main functions of FastCF are as follows: 1) approximating the motion can reduce the number of times the coupled filter is required to apply the coupled filter; 2) using envelope detection can reduce the search space of motion parameters by heuristically refining the search to increase speed; and 3) adding a fast post-processing step can compensate for possible inaccuracies introduced by the approximation. In the following sections, a conventional coupled filtering method will be briefly described, which serves as the basis for FastCF.

[0061] Coupled filtering method and "baseline method"

[0062] The coupled filtering method can, for example, be designed to compensate for feature motion decorrelation. Ultrasound image I before deformation A (X) can be modeled as:

[0063] I A (X)=Z(X)*H(X), (1)

[0064] Where Z(X) represents the scatterer, * represents convolution, and H(X) represents the point spread function of the imaging system. In addition, Z(X) and H(X) can be modeled as follows:

[0065]

[0066] here,

[0067] P={(a1, p1), (a2, p2), ..., (a N , p N )} (3)

[0068] Among them, a1,…,a N represents the amplitude of the scatterer and p1,…,p N represents the position of the scatterer, X, represent the image coordinates and the position of the scatterer respectively, N is the number of scatterers, Where x, y and z represent the horizontal direction, vertical direction and axial direction respectively, and U0=[0 0 u z ] T is the parameter of the point spread function.

[0069] Given a motion model:

[0070]

[0071] The affine motion matrix

[0072] And the displacement vector T=[t x t y t z ] T .

[0073] The following equations can be derived to describe the coupled filtering method:

[0074] I A (X)*H(MX)=I B [M(X+T)]*H(X) (5)

[0075] Among them I B (X) is the deformed image.

[0076] In particular, the coupled filtering method combines the filtering step based on (6) with block matching to compensate for feature motion decorrelation and can thus achieve high accuracy in large tissue deformation analysis. Figure 2A (right side) and Figure 2B The workflow of the coupled filtering method is shown in Figure 2. It should be noted that the motion model (5) can allow T to be searched by simply sampling across the inverse mapped image (rather than repeating the entire coupled filtering procedure). This concept can help save computational load and is adopted in the "baseline method" (described in the following paragraphs).

[0077] In addition, the coupling filtering method can adopt a coarse-to-fine search strategy to reduce the computational cost. Figure 2AAs shown in , this strategy first focuses on a small number of points and generates a coarse low-resolution (i.e., first scale / iteration) output, and then gradually refines the result by increasing the resolution and using smaller search steps in subsequent scales / iterations. Figure 2B The detailed flow of the search step is shown. The iteration is repeated until all points are processed. In this disclosure, the coupled filtering method described above and the coarse-to-fine search strategy adopted are referred to as the "baseline method".

[0078] FastCF

[0079] In this section, the present disclosure describes the algorithmic optimization of FastCF. It also describes the post-processing steps that help maintain the accuracy of the estimated strain.

[0080] FastCF parameters

[0081] First, the three search spaces are defined as follows:

[0082]

[0083] in V represents the coordinates of the voxel in the block. 范围 、V 步长 、T 范围 、 and M 范围 、 Represent the range and step size of the coordinates of the points of interest M and T respectively; Represents the element-wise less than operator; is element-wise multiplication; is a collection of vectors of size 3 by 1; is a set of vectors of size 3 by 3; is the set of integers; I3 is the 3×3 identity matrix.

[0084] In further details, Represents the point of interest. In theory, each point in the original image is associated with only one other position in the deformed image (deformed image). Therefore, only two vectors are needed to describe the two corresponding positions. Alternatively, the original position (i.e., V) and the displacement vector can also be used to describe the two corresponding positions. However, in order to determine whether the two points match, simply comparing the voxel values ​​of the two points may not be enough, because there may be other factors that may interfere with the accuracy of the results (such as noise or feature motion decorrelation). Therefore, the voxel values ​​of the voxels around the point of interest can be considered, that is, the voxel block in the original image is matched with the voxel block in the deformed image. In this regard, a simple displacement vector may not be sufficient to interpret such a large number of voxels (for example, 200 voxels in each block). Therefore, a more complex motion model is adopted, that is, an affine transformation, X′=M(X+T).

[0085] Next, this disclosure describes how to reduce the search space In order to speed up the search. From equations (1) and (6), the following equation (7) can be obtained:

[0086] I B [M(X+T)]=Z(X)*H(MX) (6)

[0087] Comparing the image before deformation described in equation (1) with the image after deformation described in equation (7), it can be seen that, among other aspects (if any), the key difference lies in H(X) and H(MX). In addition, the point spread function H(X) described in equation (3) has a Gaussian shape in both the x-direction and the y-direction, but is a Gaussian weighted cosine function (e.g., a Gabor function) in the z-direction. Therefore, by comparing H(x) with H(MX), it can be observed that changes in all three directions may cause interference between the two Gaussian envelopes, but changes only in the z-direction may cause significant interference between the peaks and troughs of the Gaussian weighted cosine function, which may lead to greater decorrelation between the image before deformation (as described in equation (1)) and the image after deformation (as described in equation (7)), and thus reduce the accuracy of the estimated hardness.

[0088] However, if you define The observation results can be summarized as follows:

[0089] When |m xx -1|+|m xy |+|m xz |+|m yx |+|m yy -1|+|m yz When |≤0.6,

[0090] H(Mc X)≈H(MX) (7)

[0091] It is worth reiterating that M c and|m xz -1|+|m xy |+|m xz |+|m yx |+|m yy -1|+|m yz |≤0.6 is formulated under the assumption that only changes in the z direction will cause significant interference between the peaks and troughs of the Gaussian weighted cosine function. In other words, as long as the above defined conditions are met, the motion parameter m xx 、m xy 、m xz 、m yx 、m yy and m yz The change in can be ignored. In addition, under the conditions defined above, it should be noted that

[0092] In further details regarding equation (8), first formulate the following equation based on both ends of equation (6):

[0093] I 1-ori (X)=I A (X)*H(MX) (9)

[0094] I 2-ori (X)=I B [M(X+T)]*H(X) (10)

[0095] However, performing convolution for all possible M based on the above two equations (9) and (10) may take too much time. Therefore, in order to reduce the amount of time required to perform convolution, the following three equations are further formulated:

[0096] I1(X)=I A (X)*H(M c X) (11)

[0097]

[0098] I2(X)=I 2t [M(X+T)] (13)

[0099] Using equations (11) and (12), I1(X) and I 2t (X) Only in M c Calculation is only required when M changes. c With only 3 independent variables (i.e., m zx、m zy and m zz ), thus the computation time can be significantly reduced compared to calculations with a maximum of 8 independent variables (9 in total, minus 1 due to the assumption of tissue incompressibility). In addition, I 2t [M(X+T)] (i.e., equation (13)). Therefore, the time complexity of calculating equation (13) is The time complexity of convolution is Therefore, it is obviously more efficient to adopt equations (11) to (13).

[0100] Furthermore, based on equations (12) and (13), the following equation (14) can be derived:

[0101]

[0102] As described above, under the definition conditions associated with equation (8), Therefore, the following relationship can be obtained:

[0103]

[0104] Similarly, based on equation (8), the following equation (16) can be derived:

[0105] when hour,

[0106] I1(X)=I A (X)*H(M c X)≈I A (X)*H(MX)=I 1-ori (X) (16)

[0107] where S is defined as the volume of the imaging region. This condition defines that the number of scatterers used in the simulation should reach approximately 10 per resolution unit to form speckles in the ultrasound image.

[0108] Since equation (6) defines I A (X)*H(MX)=I B [M(X+T)]*H(X), then I1(X)≈I2(X).

[0109] In summary, based on equations (14) to (16), the following relationships can be derived:

[0110]

[0111] In particular, as will be appreciated by those skilled in the art, the approximate equation (17) states that: A (X)*H(Mc X) and Afterwards, if the changes in the other six motion parameters satisfy the conditions defined above, i.e., are less than the thresholds defined above, there is no need to recalculate the convolution. Advantageously, since the convolution is faster than the calculation of I 2t [M(X+T)] is more time-consuming, so this implementation can help save time. In addition, the filtering scheme in different stages of the coarse-to-fine iteration can be adjusted. For example, during the last two iterations (or in other words, during the last two scales), the filtering step will account for more than 90% of the total computational load. Therefore, the user can apply the filtering step only once without affecting the accuracy of the result. In particular, since the search step size M is smaller than the search space in these scales, the filtering step size M will be smaller than the search space in these scales. 步长 is designed to be small so that such small offsets do not significantly corrupt the results.

[0112] Furthermore, in an example embodiment, three motion parameters (eg, m xy 、m yy and m zy ). This omission can be accomplished by flattening the 3D block in the y-direction (vertically), which renders these three variables / motion parameters meaningless in the calculation of M(X+T). Initially, the incompressibility constraint has omitted one motion parameter from the total number of motion parameters to be searched, i.e., reducing the total number of motion parameters from 9 to 8. However, in the example embodiment described above, the three motion parameters (i.e., m xy 、m yy and m zy ) is eliminated, which makes the incompressibility constraint no longer effective. Therefore, m xy 、m yy and m zy Omitted from a total of 9 motion parameters (instead of 8). Therefore, this embodiment can advantageously speed up the calculation by reducing the total motion parameters searched by the two motion parameters from M (i.e., 8 motion parameters (incompressibility constraint takes effect)) to 6 motion parameters. It will also be appreciated by those skilled in the art that more motion parameters can be removed depending on the desired accuracy. Additionally, in the above example embodiment, the axial direction (z direction) is retained because the axial direction has the highest resolution. In addition, the vertical direction (y direction) is omitted / flattened because the vertical direction generally has a lower resolution than the lateral direction. In other words, retaining the axial direction and the lateral direction can ensure that accurate results are obtained.

[0113] In summary, M c It can be defined as the matrix used in the filtering step as follows:

[0114]

[0115] The following Algorithm 1 presents the algorithm with the updated filtering step. The next section will describe the envelope function. In Algorithm 1, multiple loop sequences are defined for different search scales. In particular, refer to Figure 3 , in the first scale / iteration (scale / iteration=1), the search range of T (i.e. ) is relatively large, and the search space of M (i.e. ) is the same for all points of interest, so it is worth computing over the entire image. On the other hand, subsequent scales / iterations (scale / iteration > 1) have smaller And different points have different It is therefore suitable for computing only a portion of the image.

[0116]

[0117]

[0118] Envelope detection

[0119] After the filtering step as described above, the images may retain the wave shape, and their envelopes may include useful information that can be used for deformation analysis. Therefore, embodiments of the present disclosure provide methods for detecting these envelopes and aim to utilize one or more image envelopes in the filtering step so that the search can be performed at the coarsest scale with a much larger step size. In more detail, the radio frequency (RF) signal is the raw signal from the ultrasound sensor, while the B-mode signal is easier to interpret. Embodiments of the present disclosure provide different workflows to detect one or more envelopes of the filtered RF signal and convert it into one or more filtered B-mode signals, such as Figure 4 As shown in , this allows the search step size to be enlarged without causing severe peak jumping problems in RF signal-based motion tracking. In particular, referring to Figure 4 A filtering step may be applied to the RF image (e.g., the pre-deformation image or the post-deformation image). Thereafter, envelope detection may be performed on the filtered RF image. Finally, a spatial transformation step may be performed to obtain a filtered B-mode image. The order in which these processing steps are performed may be varied. Advantageously, performing the spatial transformation step after filtering the RF image may increase computational speed.

[0120] In addition, the original RF signal can have a very large spatial frequency in the axial direction. For example, for a 3 MHz ultrasound signal and a sampling rate of 15.4 MHz, peaks and troughs can be found every 5.1 voxels. This means that with a slight offset in M ​​and T, the resulting image I(MX+T) may look different and may easily mismatch the peaks and troughs. Therefore, the search steps of M and T are set relatively small so that the blocks can be correctly matched by combining the patterns of the entire block. Using envelope detection, the spatial frequency of the waveform signal can be reduced by about 5 to 7 times, as shown in Figure 2. Figure 4 As shown in , for example, there is one peak and one trough every 30 voxels. This means that the original search step size of M and T can be increased by about 6 times, and it will still enable the peaks to be correctly matched with the troughs while alleviating the peak jumping problem. Therefore, a coarse search can be applied at the beginning and the result can be finalized with a smaller search step size later. Importantly, it will be understood by those skilled in the art that the above conversion process can be applied together with a coarse-to-fine multi-scale strategy to further accelerate it. The envelope is the L1 norm of the Hilbert transform signal. The updated algorithm for the coarse-to-fine strategy is shown in Algorithm 2.

[0121] It is worth mentioning that flattening is performed on blocks (or one or more voxels) of the image rather than on the image itself. In other words, when using a similarity metric to measure the similarity between two corresponding points in the pre-deformed image and the deformed image, only the flattened region of interest (i.e., a block containing X*1*Z voxels) is considered in both the pre-deformed image and the deformed image. In particular, when calculating the similarity metric, only the voxels in the 2D block 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 in these 2D blocks, because voxels in other vertical planes may affect the accuracy of the results.

[0122]

[0123] Post-processing

[0124] Peak jumps are common in tissue deformation analysis based on RF images and can degrade motion tracking results. To address this issue, embodiments of the present disclosure provide a post-processing step based on majority voting. In particular, it can be observed that the axial displacement (i.e., D z ) is more accurate than the other two directions. Therefore, based on this observation, D can be calculated in three directions. z For example, by using PQ = [V step，x 0 0] T Define two adjacent points p, And use the Lagrangian motion model to define the displacement D opt (V)=Mopt (V)(V+T opt (V))-V, the following restrictions can be derived:

[0125]

[0126] Where [m zx，min , m zx，max ] is m zx Search scope.

[0127] Similar restrictions can be established along the other two directions. Figure 5 , any gradient exceeding these theoretical limits may be considered erroneous, and voxels surrounded by erroneous gradients may be defined as erroneous voxels. Subsequently, the M and T of the erroneous voxels can be interpolated based on the correct voxels, and the results can be refined at a finer scale. Typically less than 5% are erroneous voxels, and the other voxels are unaffected. In addition, the post-processing steps can be designed to be lightweight and not significantly slow down the overall computation time. The results of the post-processing steps will be discussed in later sections of this disclosure.

[0128] The post-processing step algorithm is shown in Algorithm 3 below.

[0129]

[0130] FastCF algorithm

[0131] The overall algorithm of FastCF is presented in Algorithm 4 and described in this section. When searching for the optimal T in Algorithm 1, a summation table approach can be adopted. In the example embodiment, the CPU implementation is coded in MATLAB. In the first scale / iteration, after calculating I1(X) and I 2temp After (X), the calculation is parallelized; in the remaining scales, all V calculations are parallelized. In addition, the MATLAB code for interpolation and optimal value search using C++MEX functions is reimplemented. Pointers and dimension permutations are also used. In the post-processing step, the graph can be stored in the form of a 3D matrix, where elements with odd coordinates represent vertices, such as Figure 5 Each vertex represents A voxel in , and only if their difference D opt Only when it falls within the theoretical limit can it be connected to its adjacent vertices.Then the MATLAB built-in function bwlabeln(), which helps to find connected components in binary images, is used to complete the relevant graph operations.

[0132]

[0133] Figure 6An overview of an exemplary process flow of the method described in this disclosure is provided. In other words, the embodiments described in this disclosure provide an accelerated coupled filtering method for tissue deformation analysis. This method can be used Figure 7 The system 700 shown is implemented as follows, Figure 7 A schematic diagram of a system 700 for performing an accelerated coupled filtering method for tissue deformation analysis is shown. The system 700 may include a processing device 702. In an embodiment of the present disclosure, the system 700 may be communicatively coupled to an ultrasound system 704 to receive and transmit information / data (e.g., ultrasound images, etc.). However, it will be appreciated by those skilled in the art that, depending on the application, the system 700 may be communicatively coupled to other suitable systems or devices for performing the method.

[0134] refer to Figure 8 , the method 800 may include the following steps:

[0135] Step 802: Applying a first filter to the pre-deformation image of the tissue by a processing device to obtain a filtered pre-deformation image.

[0136] Step 804: Applying, by the processing device, a second filter to the deformed image of the tissue to obtain a filtered deformed image.

[0137] The filtered pre-deformation image and the filtered post-deformation image are correlated using a first motion matrix comprising a plurality of first motion parameters. The plurality of first motion parameters comprises at least three basic first motion parameters. Each of the three basic first motion parameters represents movement of the tissue along the axial direction relative to the axial direction, the vertical direction, and the transverse direction, respectively.

[0138] Step 806: Estimate, by the processing device, a corresponding value of each of the plurality of first motion parameters. Each of the estimated corresponding values ​​represents a difference between the filtered image before deformation and the filtered image after deformation.

[0139] Step 808: In response to determining that at least one of the estimated corresponding values ​​satisfies at least one predefined criterion, updating at least one of the estimated corresponding values ​​to an optimal value of at least one corresponding first motion parameter among the plurality of first motion parameters. Those skilled in the art will appreciate that the term "corresponding" defines, by way of example, that if only the estimated corresponding values ​​of a primary motion parameter representing movement of tissue along the axial direction relative to the transverse direction satisfy at least one predefined criterion, then only the corresponding values ​​of such primary motion parameters will be updated. The remaining primary motion parameters will not be updated.

[0140] In an embodiment of the present disclosure, method 800 may include the following steps: detecting, by a processing device, one or more envelopes present in a filtered pre-deformation image and a filtered post-deformation image to obtain a filtered B-mode pre-deformation image and a filtered B-mode post-deformation image.

[0141] In an embodiment of the present disclosure, applying steps 802 and 804 may include convolving a first point spread function and a second point spread function of the ultrasound system with the pre-deformed image and the deformed image, respectively. The first point spread function may be a modified version of the second point spread function.

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

[0143] In embodiments of the present disclosure, the filtered pre-deformation image and the filtered post-deformation image may be three-dimensional (3D) images. Those skilled in the art will appreciate that the pre-deformation image, the post-deformation image, the filtered pre-deformation image, and the filtered post-deformation image are not limited to 3D images and may encompass other image formats (such as two-dimensional images) depending on user preferences and / or applications. Furthermore, estimating the corresponding value may include the following sub-steps.

[0144] Sub-step A: Flattening, by a processing device, a first plurality of voxels and a second plurality of voxels around each of the one or more pre-deformation points of interest in the filtered pre-deformation image and around each of the one or more post-deformation points of interest in the filtered post-deformation image. The flattening may be performed along a vertical direction of each of the first and second pluralities of voxels, which may reduce the first and second motion parameters.

[0145] In an embodiment of the present disclosure, the estimation step 806 can be performed in multiple iterations. In the first iteration of the multiple iterations, the multiple first motion parameters can include two basic first motion parameters of the three basic first motion parameters that respectively represent the movement of the tissue along the axial direction relative to the axial direction and the transverse direction, and the remaining basic first motion parameters can be removed by flattening. Between the second iteration and the final iteration of the multiple iterations, the multiple first motion parameters can include two basic first motion parameters of the three basic first motion parameters, two transverse first motion parameters that respectively represent the movement of the tissue along the transverse direction relative to the transverse direction and the axial direction, and two vertical first motion parameters that respectively represent the movement of the tissue along the vertical direction relative to the transverse direction and the axial direction.

[0146] In an embodiment of the present disclosure, the second motion matrix may include a plurality of second motion parameters, and the estimation step 806 may further include the following sub-steps.

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

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

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

[0150] Sub-step E: generating, by the processing device, one or more deformed voxels based on the search value of each of the plurality of second motion parameters.

[0151] Sub-step F: Calculating, by the processing device, a similarity index based on one or more first voxels surrounding the one or more pre-deformed voxels and one or more second voxels surrounding the one or more post-deformed voxels. The similarity index may define a similarity between the target pre-deformed interest point and the post-deformed interest point corresponding to the target pre-deformed interest point.

[0152] Sub-step G: The processing device determines whether the target interest point before deformation matches the corresponding interest point after deformation based on the calculated similarity index.

[0153] Additionally, in an embodiment of the present disclosure, the plurality of second motion parameters may include at least two basic second motion parameters, two lateral second motion parameters, and two vertical motion parameters. The estimation step 806 may further include, in a first iteration, determining, by the processing device, a value for each of two basic first motion parameters among the three basic first motion parameters based on steps (A) to (G). Between the second iteration and the third final iteration, determining, by the processing device, an updated value for each of two basic first motion parameters among the three basic first motion parameters based on steps (A) to (G), and assigning, by the processing device, the value of each of the two lateral second motion parameters and the two vertical second motion parameters searched for in step (C) in the previous iteration as the corresponding value of each of the two lateral first motion parameters and the two vertical first motion parameters. In the final two iterations, the determined updated value of each of the two basic first motion parameters of the three basic first motion parameters in the third final iteration is assigned by the processing device as the optimal value of each of the two basic first motion parameters of the three basic first motion parameters, and the updated value of each of the two horizontal second motion parameters and the two vertical second motion parameters searched in step (C) in the third final iteration is assigned by the processing device as the optimal value of each of the two horizontal first motion parameters and the two vertical first motion parameters.

[0154] In an embodiment of the present disclosure, the similarity index may include any of the following: normalized correlation coefficient, sum of absolute differences, and sum of squared differences. However, those skilled in the art will readily appreciate that other suitable indicators may be used depending on user preferences and / or applications.

[0155] Furthermore, method 800 may include the following steps: applying, by the processing device, upper and lower limits to the gradient of the displacement range of the fundamental first motion parameter along the axial direction. Values ​​falling outside the displacement range may be determined as erroneous values ​​of the fundamental first motion parameter. In response to determining that one or more voxels include erroneous values ​​of the fundamental first motion parameter, the processing device may determine, by the processing device, corrected values ​​of the fundamental first motion parameter for the one or more voxels. The corrected values ​​may be determined by interpolating the erroneous values ​​based on predetermined corrected values ​​of the fundamental first motion parameter for voxels surrounding the one or more voxels.

[0156] Evaluate

[0157] Experimental setup

[0158] In this section, we evaluate the timing performance and tracking accuracy of FastCF. The experimental setup includes hardware, software, two simulated datasets and one real dataset, a set of parameters, and an evaluation metric. The CPU platform is an AMD Ryzen Threadripper PRO 3995WX (64 cores) with 512GB of memory. The CPU program utilizes less than 64GB of memory with 64 MATLAB worker threads. The CPU code was implemented and tested in MATLAB R2021b on Windows Server 2022.

[0159] Dataset

[0160] To simulate a dataset with known ground truth, elastic tissue was modeled using SolidWorks. Figure 9 As shown in , the elastic tissue has hard spherical regions fused in a soft cubic material. The model was then uniformly compressed by 5% in the axial direction, and the movement under compression was simulated using SolidWorks. In addition, based on the simulation results, a pair of 3D ultrasound images were generated using the simulation software Field II, and the parameters of the point spread function H(X) were estimated. The ultrasound images were stored in RF format, and the image size was 101x101x1001 voxels. Additionally, a real 3D dataset was collected using a Philips iE33 xMATRIX ultrasound system. 3D B-mode images of the subject's left ventricle were acquired using an X3-1 transducer, and the image size was 302×158×302 voxels. Two volumes were then manually selected during contraction for illustration. Since the RF signal is not accessible, an approximate version of the coupled filtering step is adopted by removing the filtering step, as shown below:

[0161] I A (X)≈I B [M(X+T)] (20)

[0162] In addition, FastCF is tested using a 2D simulated dataset and the results are compared with baseline methods.

[0163] Parameter settings

[0164] Figure 10 The parameters and values ​​of the simulated 3D dataset are defined in . In addition, Figure 10 The time complexity, space complexity, and estimated floating point operations per second (FLOPs) for 2D and 3D datasets are also shown. Figure 10 In , (a) represents coupled filtering with exhaustive search, (b) represents a coarse-to-fine search strategy, (c) represents an approximation of the matrix M, (d) represents envelope detection, and (e) represents a summation table. Other parameters are in Figure 3 In addition, Figure 11 The parameters for time complexity and space complexity analysis of 2D and 3D implementations of various coupled filtering methods are shown. xx and m zz The search range is between -10% and 10%, and m xz 、m yx 、m yz 、m zx The search range for is between -5% and 5%. These search ranges are large enough for the dataset.

[0165] Evaluation Metrics

[0166] In the experiment, the square root of the mean square error (SRMSE) of the axial strain is used to measure the accuracy of the estimation results:

[0167]

[0168] where ∈(W)=m zz -1 is the estimated axial strain of a certain block B, ∈0(B) represents the corresponding ground truth, and N B is the number of blocks.

[0169] In the MATLAB program, the timer "tic toc" was used to measure the running time, and "mpiprofile" was used to measure the running time of different functions. Only the time to initialize the MATLAB parallel pool was excluded.

[0170] FLOP estimation and complexity analysis

[0171] Figure 12 The detailed time complexity of all relevant methods is shown. As an example, this disclosure also demonstrates how to calculate the time complexity and FLOPs of exhaustive search.

[0172] For the filtering step, the size is The kernel and size are The image is convolved, which can be accelerated by the circular convolution property of Fourier transform. In particular, a fast Fourier transform (FFT) is applied to the padded image, and the result is multiplied by the transformed kernel, and then an inverse FFT is applied, and the final result is obtained, that is, 4 forward FFTs and 2 inverse FFTs to obtain two convolution results I A (X)*H(MX) and I B [M(X+T)]*H(X). Since I A(X) and H(X) remain unchanged, so the entire procedure requires only 2 forward FFTs. Meanwhile, a 1D FFT with data size 1x n requires 5n log2n operations. 3D FFT shares a similar definition except for the three nested summations, and each summation has the same complexity as the 1D FFT, so a 3D FFT requires 15n 3 log2n 3 The inverse FFT also has the same operations as the FFT, and the multiplication is performed n times. 3 operations. Therefore, the total FLOPs for the filtering step is [15(n i +n p ) 3 log(n i +n p ) 3 +2(n i +n p ) 3 ], and the time complexity is

[0173] Furthermore, a block matching algorithm is employed to search for the best result. For the interpolation step, interpolation is performed only on each M image, and T can be implemented by sliding a window in the interpolated image. This implementation requires 30 operations per voxel, so the total FLOPs is And the time complexity is As for the indicator calculation, for the output Voxel, through and The dimensions around it are The search is performed on blocks of , and 3 additions and 3 multiplications are performed on each voxel. There are also some operations per index, but they are negligible compared to the other indexes. Therefore, the rough FLOPs for block matching is And the time complexity is In summary, the total time complexity is

[0174] Furthermore, by taking the approximate value of M (as described above) and referring to Figure 12 , the convolution can be obtained from times reduced to And used to calculate the time I of affine motion 2t [M(X+T)] can be obtained from Reduce to Therefore, the time complexity can be reduced to If these n a blocks, the space complexity can be turned into Furthermore, if envelope detection is applied, the approximation of the matrix M can reduce the number of times the affine motion is calculated from is reduced to a constant, so the main computation load is offset convolution, and the time complexity is reduced to

[0175] Results on 3D simulated datasets

[0176] On the simulated 3D dataset, Figure 10 The estimated FLOPs and complexity of different methods are shown in Figure 2. Both the M-approximation and envelope detection methods effectively reduce the overhead, and the use of the summation table method also helps. In summary, FastCF is estimated to reduce 4.21 billion FLOPs compared to the baseline method. Figure 13 The average SRMSE and computation time for each scale are shown in . Figure 14 The results are shown in . Figure 13 As shown in , as the search step size of M and T becomes smaller, the SRMSE decreases continuously. Figure 14 , the final result is in the cross-section diagram ( Figure 14 (a) and Figure 14 (c)) clearly shows the hard sphere, achieving the goal of tissue deformation analysis. When carefully observing the running times at different scales, it is observed that, with the exception of the first scale (which initializes the search), each of the other scales shows an increase in the number of blocks by approximately eight times. Therefore, the gradual reduction in the search space is not enough to compensate for the significant increase in the number of blocks, which leads to an increase in the overall computational time. Additionally, as described above, the filtering strategy for the fifth scale is adjusted. Therefore, the fifth scale takes a shorter time than the other scales. Because the interpolation between the scales takes a few seconds, the total time is slightly longer than the sum of all six scale times.

[0177] Furthermore, the computation time of different functions is analyzed, e.g. Figure 15 In this analysis, the program is divided into four main components, and the first two main components are further decomposed. Figure 15 As shown in Figure 2, the coupled filter component is the most time-consuming, with the Fourier transform relying on built-in libraries and unable to be optimized externally. Data preparation and frequency domain multiplication may be limited by memory speed. In addition, the amount of time spent on the post-processing component is negligible compared to the search process. In the CPU implementation, the long computation time of other functions is mostly spent on MATLAB indexing.

[0178] Results on 3D real-world datasets

[0179] For the real 3D dataset, the white area near the ventricular wall was manually segmented in advance. Figure 16This analysis required approximately 295 seconds in FastCF to generate a strain image of 97 × 65 × 129 pixels on a CPU. Figure 16 , the overall compression trend in the slices is clearly visible, and further interpretation of these results will require clinical information. This dataset was generated at a frame rate of only 8 Hz. In other words, embodiments of the present disclosure can produce good results with ultrasound images sampled at a lower frame rate (e.g., 8 Hz). In contrast, traditional deformation analysis methods require a frame rate of at least 40 Hz to 80 Hz to produce good results.

[0180] Comparison with Baseline Methods

[0181] For comparison, FastCF and the baseline method are performed on 2D simulated data including a hard cylinder. This is because executing the baseline method on 3D data would take too much time. Figure 14 (a) Similar. Figure 17 As shown in Figure 2, both the baseline method and FastCF can identify hard inclusion compounds in a 2D dataset. Furthermore, FastCF achieves an SRMSE of 1.11%, while the baseline method achieves an SRMSE of 1.40%. FastCF also completes the analysis in approximately 12.9 seconds on a CPU, while the baseline method takes approximately 1.72 days on the same CPU and 20 minutes on an FPGA.

[0182] Verification of the approximate value of the matrix M

[0183] In this subsection, the effects of three parameters in M ​​(defined in Equation (5)) are studied. First, the image model in Equation (1) is used and the scatter density requirement defined in Equation (8) is applied to generate pre- and post-deformation images with different motion matrices M. Then, 27 blocks are uniformly selected and their normalized correlation coefficients (NCCs) are averaged to plot the NCCs against the nine motion parameters in M, respectively. The results are given in Figure 18 Reference Figure 18 , it is observed that the result is consistent with equation (8), that is, if along m x* and m y* There is only a small offset, and m z* unchanged, then I A (X) and I B [M(X+T)] can still achieve a high correlation coefficient.

[0184] In addition, in this experiment, four options (i.e., "very slow" to "medium") were tested without post-processing to study the effect of different search spaces for M and T. The runtimes for these options were 434.607s, 83.928s, 15.796s, 2.571s, and 2.589s, respectively. For comparison purposes, FastCF, i.e., the "medium" setting with post-processing, was also included in this experiment. Figure 19 The errors across six scales are shown using the five options. Compared to the "Medium" setting with no post-processing, the slower options achieve lower SRMSE but show diminishing returns. On the other hand, FastCF is able to achieve similar accuracy to the other slower options and runs almost as fast as the "Medium" option.

[0185] Figure 20 Depicted is an exemplary computing device 2000, hereinafter interchangeably referred to as computer system 2000. The following description of computing device 2000 is provided by way of example only and is not intended to be limiting.

[0186] like Figure 20 As shown, the example computing device 2000 includes a processor 2002 for executing software routines. Although a single processor is shown for clarity, the computing device 2000 may also include a multi-processor system. The processor 2002 is connected to a communication infrastructure 2004 for communicating with other components of the computing device 2000. The communication infrastructure 2004 may include, for example, a communication bus, a crossbar switch, or a network.

[0187] The computing device 2000 also includes a main memory 2006 (such as random access memory (RAM)) and a secondary memory 2008. The secondary memory 2008 may include, for example, a hard disk drive 2010 and / or a removable storage drive 2012, which may include a floppy disk drive, a magnetic tape drive, an optical disk drive, or the like. The removable storage drive 2012 reads from and / or writes to a removable storage unit 2014 in a well-known manner. The removable storage unit 2014 may include a floppy disk, a magnetic tape, an optical disk, or the like, which is read from and written to by the removable storage drive 2012. As will be appreciated by those skilled in the relevant art, the removable storage unit 2014 includes a computer-readable storage medium having stored therein computer-executable program code instructions and / or data.

[0188] In alternative embodiments, the secondary memory 2008 may additionally or alternatively include other similar means for allowing computer programs or other instructions to be loaded into the computing device 2000. Such means may include, for example, a removable storage unit 2016 and an interface 2018. Examples of removable storage units 2016 and interfaces 2018 include program cartridges and cartridge interfaces (such as those found in video game console devices) that allow software and data to be transferred from the removable storage unit 2016 to the computer system 2000, removable storage chips (such as EPROM or PROM) and associated sockets, and other removable storage units 2016 and interfaces 2018.

[0189] The computing device 2000 also includes at least one communication interface 2020. The communication interface 2020 allows software and data to be transferred between the computing device 2000 and external devices via a communication path 2022. In various embodiments of the present invention, the communication interface 2020 allows data to be transferred between the computing device 2000 and a data communication network (such as a public or private data communication network). The communication interface 2020 can be used to exchange data between different computing devices 2000, such computing devices 2000 forming part of an interconnected computer network. Examples of the communication interface 2020 may 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 2020 may be wired or wireless. The software and data transferred via the communication interface 2020 are in the form of signals, which may be electronic, electromagnetic, optical, or other signals that can be received by the communication interface 2020. These signals are provided to the communication interface via the communication path 2022.

[0190] like Figure 20 As shown, computing device 2000 also includes a display interface 2024 that performs operations for rendering images to an associated display 2026 and an audio interface 2028 that performs operations for playing audio content via associated speakers 2030 .

[0191] As used herein, the term "computer program product" may refer in part to a removable storage unit 2014, a removable storage unit 2016, a hard disk installed in the hard disk drive 2010, or a carrier that carries the software via a communication path 2022 (wireless link or cable) to the communication interface 2020. A computer-readable storage medium is any non-transitory tangible storage medium that provides recorded instructions and / or data to the computing device 2000 for execution and / or processing. Examples of such storage media include floppy disks, magnetic tapes, CD-ROMs, DVDs, Blu-ray discs, and the like. TMA disk, a hard drive, a ROM or integrated circuit, a USB memory, a magneto-optical disk, or a computer-readable card (such as a PCMCIA card), etc., whether such a device is internal or external to the computing device 2000. Examples of transitory or non-tangible computer-readable transmission media that may also participate in providing software, applications, instructions and / or data to the computing device 2000 include radio or infrared transmission channels and a network connection to another computer or networked device, as well as the Internet or an intranet including electronic mail transmissions and information recorded on websites, etc.

[0192] A computer program (also referred to as computer program code) is stored in the main memory 2006 and / or the secondary memory 2008. The computer program may also be received via the communication interface 2020. Such a computer program, when executed, enables the computing device 2000 to perform one or more features of the embodiments discussed herein. In various embodiments, the computer program, when executed, enables the processor 2002 to perform the features of the various embodiments described above. Thus, such a computer program represents a controller of the computer system 2000.

[0193] The software may be stored in a computer program product and loaded into the computing device 2000 using the removable storage drive 2012, the hard drive 2010, or the interface 2018. Alternatively, the computer program product may be downloaded to the computer system 2000 via the communication path 2022. The software, when executed by the processor 2002, causes the computing device 2000 to perform the functions of the embodiments described herein.

[0194] It should be understood that Figure 20 The embodiments are presented by way of example only. Thus, in some embodiments, one or more features of computing device 2000 may be omitted. Furthermore, in some embodiments, one or more features of computing device 2000 may be combined. Additionally, in some embodiments, one or more features of computing device 2000 may be separated into one or more component parts.

[0195] It is understandable that Figure 20 The elements illustrated in the figure serve to provide means for executing various functions and operations of the server as described in the above embodiments.

[0196] In one embodiment, a server can be generally described as a physical device including at least one processor and at least one memory containing computer program code. The at least one memory and the computer program code are configured to work with the at least one processor to cause the physical device to perform necessary operations.

[0197] When the computing device 2000 is configured to implement the system 700 for performing the acceleration coupled filtering method for tissue deformation analysis, the system 700 may have a non-transitory computer-readable medium having stored thereon an application program that, when executed, causes the system 700 to perform the acceleration coupled filtering method 800 described above and / or any other method described herein.

[0198] It will be appreciated by those skilled in the art that, without departing from the spirit or scope of the invention as broadly described, various changes and / or modifications may be made to the invention shown in the specific embodiments. Therefore, the present embodiments are to be considered in all respects as illustrative and not restrictive.

Claims

1. An accelerated coupled filtering method for tissue deformation analysis, comprising: applying, by a processing device, a first filter on the pre-deformation image of the tissue to obtain a filtered pre-deformation image; applying, by the processing device, a second filter on the deformed image of the tissue to obtain a filtered deformed image, wherein the filtered pre-deformed image and the filtered deformed image are related by a first motion matrix comprising a plurality of first motion parameters, and wherein the plurality of first motion parameters comprises at least three basic first motion parameters, each of the three basic first motion parameters representing movement of the tissue along an axial direction relative to an axial direction, a vertical direction, and a lateral direction, respectively; estimating, by the 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 filtered image before the deformation and the filtered image after the deformation; as well as In response to determining that at least one of the estimated corresponding values ​​satisfies at least one predefined criterion, at least one of the estimated corresponding values ​​is updated to an optimal value of at least one corresponding first motion parameter of the plurality of first motion parameters.

2. The method according to claim 1, further comprising: One or more envelopes present in the filtered pre-warping image and the filtered post-warping image are detected by the processing device to obtain a filtered B-mode pre-warping image and a filtered B-mode post-warping image.

3. The method of claim 1 , wherein applying the first filter and the second filter comprises convolving a first point spread function and a second point spread function of the 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.

4. The method of claim 1 , wherein the method further comprises, after applying the second filter, spatially transforming, by the processing device, the warped image based on at least a second motion matrix, and wherein the first motion matrix is ​​a modified version of the second motion matrix.

5. The method of claim 4 , wherein the filtered pre-deformation image and the filtered post-deformation image are three-dimensional (3D) images, and wherein estimating the corresponding value further comprises: (A) flattening, by the processing device, a first plurality of voxels and a second plurality of voxels, the first plurality of voxels and the second plurality of voxels respectively around each of one or more pre-deformation points of interest in the filtered pre-deformation image and around each of one or more post-deformation points of interest in the filtered post-deformation image, wherein the flattening is performed along a vertical direction of each of the first plurality of voxels and the second plurality of voxels, and wherein the flattening reduces the plurality of first motion parameters and the plurality of second motion parameters.

6. The method of claim 5 , wherein estimating the corresponding value is performed in a plurality of iterations, and wherein: In a first iteration of the plurality of iterations: The plurality of first motion parameters include two of the three basic first motion parameters representing movement of the tissue along the axial direction relative to the axial direction and the transverse direction, respectively, and wherein the remaining basic first motion parameters are removed by flattening, and Between the second iteration and the final iteration of the plurality of iterations: The multiple first motion parameters include two of the three basic first motion parameters, the two lateral first motion parameters respectively representing the movement of the tissue along the lateral direction relative to the lateral direction and the axial direction, and the two vertical first motion parameters respectively representing the movement of the tissue along the vertical direction relative to the lateral direction and the axial direction.

7. The method of claim 6, wherein the second motion matrix comprises a plurality of second motion parameters, and wherein estimating the corresponding values ​​further comprises: (B) searching, by the processing device, for each of the one or more pre-deformation interest points in the filtered pre-deformation image, for a value of each of the plurality of first motion parameters; (C) searching, by the processing device, for each of the one or more warped interest points in the filtered warped image, for a value of each of the plurality of second motion parameters; (D) generating, by the processing device, one or more pre-deformation voxels based on the search value of each of the plurality of first motion parameters; (E) generating, by the processing device, one or more deformed voxels based on the search value of each of the plurality of second motion parameters; (F) calculating, by the processing device, a similarity index 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 index defines a similarity between a target pre-deformation interest point and a post-deformation interest point corresponding to the target pre-deformation interest point; as well as (G) Determining, by the processing device based on the calculated similarity index, whether the target pre-deformation interest point matches the corresponding post-deformation interest point.

8. The method of claim 7 , wherein the plurality of second motion parameters comprises at least two basic second motion parameters, two lateral second motion parameters, and two vertical second motion parameters, and wherein estimating the corresponding values ​​comprises: In the first iteration: determining, by the processing device, a value for each of two of the three basic first motion parameters based on steps (A) to (G); Between the second iteration and the third final iteration: determining, by the processing device, updated values ​​for each of two of the three basic first motion parameters based on steps (A) to (G); and assigning, by the processing device, a value of each of the two lateral second motion parameters and the two vertical second motion parameters searched for in step (C) in a previous iteration as a corresponding value of each of the two lateral first motion parameters and the two vertical first motion parameters; as well as In the final two iterations: assigning, by the processing device, the updated value for each of the two of the three basic first motion parameters determined in the third final iteration as the optimal value for each of the two of the three basic first motion parameters; and The updated value of each of the two lateral second motion parameters and the two vertical second motion parameters searched in step (C) in the third final iteration is assigned by the processing device as the optimal value of each of the two lateral first motion parameters and the two vertical first motion parameters. 9 . The method according to claim 8 , wherein the similarity index comprises any one of a normalized correlation coefficient, a sum of absolute differences, and a sum of squared differences.

10. The method according to claim 1, further comprising: imposing, by the processing device, upper and lower limits on a gradient of a displacement range of the basic first motion parameter along the axial direction, wherein values ​​falling outside the displacement range are determined to be erroneous values ​​of the basic first motion parameter; In response to determining that the one or more voxels include an erroneous value of the base first motion parameter: A correction value for the base first motion parameter of the one or more voxels is determined by the processing device, wherein the correction value is determined by interpolating the error value based on predetermined correction values ​​for the base first motion parameter of voxels surrounding the one or more voxels.

11. A system for performing an acceleration coupled filtering method for tissue deformation analysis, the system comprising a processing device configured to: applying a first filter to the pre-deformation image of the tissue to obtain a filtered pre-deformation image; applying a second filter to the deformed image of the tissue to obtain a filtered deformed image, wherein the filtered pre-deformed image and the filtered deformed image are correlated by a first motion matrix comprising a plurality of first motion parameters, and wherein the plurality of first motion parameters comprises at least three basic first motion parameters, each of the three basic first motion parameters representing movement of the tissue along an axial direction relative to an axial direction, a vertical direction, and a lateral direction, respectively; estimating 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-warping filtered image and the post-warping filtered image; as well as In response to determining that at least one of the estimated corresponding values ​​satisfies at least one predefined criterion, at least one of the estimated corresponding values ​​is updated to an optimal value of at least one corresponding first motion parameter of the plurality of first motion parameters.

12. The system of claim 11, further comprising: One or more envelopes present in the filtered pre-deformation image and the filtered post-deformation image are detected to obtain a filtered B-mode pre-deformation image and a filtered B-mode post-deformation image.

13. The system of claim 11 , wherein to apply the first filter and the second filter, the processing device is configured to convolve a first point spread function and a 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.

14. The system of claim 11, wherein the processing device is further configured to, after applying the second filter, spatially transform the warped image based on at least a second motion matrix, and wherein the first motion matrix is ​​a modified version of the second motion matrix.

15. The system of claim 14 , wherein the filtered pre-deformation image and the filtered post-deformation image are three-dimensional (3D) images, and wherein to estimate the corresponding value, the processing device is configured to: (A) flattening a first plurality of voxels and a second plurality of voxels, the first plurality of voxels and the second plurality of voxels respectively being around each of one or more pre-deformation points of interest in the filtered pre-deformation image and around each of one or more post-deformation points of interest in the filtered post-deformation image, wherein the flattening is performed along a vertical direction of each of the first plurality of voxels and the second plurality of voxels, and wherein the flattening reduces the plurality of first motion parameters and the plurality of second motion parameters.

16. The system of claim 15, wherein estimating the corresponding value is performed in a plurality of iterations, and wherein: In a first iteration of the plurality of iterations: The plurality of first motion parameters include two basic first motion parameters of the three basic first motion parameters respectively representing movement of the tissue along the axial direction relative to the axial direction and the transverse direction, and wherein the remaining basic first motion parameters are removed by the flattening, and Between the second iteration and the final iteration of the plurality of iterations: The multiple first motion parameters include the two basic first motion parameters of the three basic first motion parameters, two lateral first motion parameters respectively representing the movement of the tissue along the lateral direction relative to the lateral direction and the axial direction, and two vertical first motion parameters respectively representing the movement of the tissue along the vertical direction relative to the lateral direction and the axial direction.

17. The system of claim 16, wherein the second motion matrix comprises a plurality of second motion parameters, and wherein to estimate the corresponding values, the processing device is further configured to: (B) searching for a value of each of the plurality of first motion parameters for each of the one or more pre-deformation interest points in the filtered pre-deformation image; (C) searching for a value of each of the plurality of second motion parameters for each of the one or more warped interest points in the filtered warped image; (D) generating one or more pre-deformation voxels based on the search value of each of the plurality of first motion parameters; (E) generating one or more deformed voxels based on the search value of each of the plurality of second motion parameters; (F) calculating a similarity index 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 index defines a similarity between a target pre-deformation interest point and a post-deformation interest point corresponding to the target pre-deformation interest point; as well as (G) Determining whether the target pre-deformation interest point matches the corresponding post-deformation interest point based on the calculated similarity index.

18. The system of claim 17 , wherein the plurality of second motion parameters comprises at least two basic second motion parameters, two lateral second motion parameters, and two vertical motion parameters, and wherein in order to estimate the corresponding values, the processing device is further configured to: In the first iteration: determining a value for each of the two basic first motion parameters of the three basic first motion parameters based on steps (A) to (G); Between the second iteration and the third final iteration: determining an updated value for each of the two basic first motion parameters of the three basic first motion parameters based on steps (A) to (G); and assigning a value of each of the two lateral second motion parameters and the two vertical second motion parameters searched for in step (C) in a previous iteration as a corresponding value of each of the two lateral first motion parameters and the two vertical first motion parameters; as well as In the final two iterations: assigning the determined updated value of each of the two of the three basic first motion parameters in a third final iteration as an optimal value for each of the two of the three basic first motion parameters; and The updated value of each of the two lateral second motion parameters and the two vertical second motion parameters searched in step (C) in the third final iteration is assigned as the optimal value of each of the two lateral first motion parameters and the two vertical first motion parameters.

19. The system of claim 18, wherein the similarity index comprises any one of a normalized correlation coefficient, a sum of absolute differences, and a sum of squared differences.

20. The system of claim 11, wherein the processing device is further configured to: imposing upper and lower limits on a gradient of a displacement range of the basic first motion parameter along the axial direction, wherein values ​​falling outside the displacement range are determined to be erroneous values ​​of the basic first motion parameter; In response to determining that one or more voxels include an erroneous value of the base first motion parameter: A correction value of the base first motion parameter of the one or more voxels is determined, wherein the correction value is determined by interpolating the erroneous value based on predetermined correct values ​​of the base first motion parameter of voxels surrounding the one or more voxels.

Images