WO2015032064A1 - 声辐射力脉冲成像估算方法和系统、计算机存储介质 - Google Patents

声辐射力脉冲成像估算方法和系统、计算机存储介质 Download PDF

Info

Publication number
WO2015032064A1
WO2015032064A1 PCT/CN2013/083049 CN2013083049W WO2015032064A1 WO 2015032064 A1 WO2015032064 A1 WO 2015032064A1 CN 2013083049 W CN2013083049 W CN 2013083049W WO 2015032064 A1 WO2015032064 A1 WO 2015032064A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
data
displacement
filtering
acoustic radiation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/CN2013/083049
Other languages
English (en)
French (fr)
Inventor
郑海荣
冯歌
王丛知
曾成志
杨戈
曾博
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenzhen Institute of Advanced Technology of CAS
Original Assignee
Shenzhen Institute of Advanced Technology of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shenzhen Institute of Advanced Technology of CAS filed Critical Shenzhen Institute of Advanced Technology of CAS
Priority to PCT/CN2013/083049 priority Critical patent/WO2015032064A1/zh
Publication of WO2015032064A1 publication Critical patent/WO2015032064A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/48—Diagnostic techniques
    • A61B8/485—Diagnostic techniques involving measuring strain or elastic properties
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52023—Details of receivers
    • G01S7/52036—Details of receivers using analysis of echo signal for target characterisation
    • G01S7/52042—Details of receivers using analysis of echo signal for target characterisation determining elastic properties of the propagation medium or of the reflective target

Definitions

  • the present invention relates to the field of medical imaging, and in particular to a method and system for estimating acoustic radiation force pulse imaging, and a computer storage medium.
  • ARFI Acoustic Radiation Force Impulse
  • Acoustic radiation force pulse imaging calculation refers to estimating the displacement of shear waves inside the tissue, and later understands the organization by judging the displacement.
  • There are various methods for estimating displacement such as real-time ultrasound imaging monitoring methods based on general-purpose GPU (Graphics Process Unit) programming.
  • the steps include: 1) expanding data; 2) windowing data; 3) cross-correlation processing data; Estimating the displacement; 5) repeating steps 3 and 4 over the entire depth of time; 6) median filtering; 7) estimating the rate.
  • the monitoring method uses a common time-domain cross-correlation algorithm, and the estimation sensitivity for small displacement is low and the accuracy is low.
  • An acoustic radiation force pulse imaging estimation method comprising:
  • the data filtering transformation step performs forward filtering processing and inverse filtering processing on each data matrix, and performs Hilbert transform on the filtered data matrix to obtain analytical data;
  • the displacement combining step performs time accumulation on each displacement rectangle, takes the cumulative and maximum column displacement, and combines the column displacements extracted by the predetermined number of displacement matrices to generate a new displacement matrix;
  • Interpolation step performing cubic spline interpolation in the time direction on the obtained new displacement matrix; fitting step, performing Radon transformation on the interpolation result to obtain a Radon transformation matrix;
  • the estimating step acquires parameters of the acoustic radiation force pulse imaging, and obtains a shear wave rate according to the parameter and the Ladong transform matrix.
  • An acoustic radiation force pulse imaging estimation system comprising:
  • a data dividing module configured to read input data collected by acoustic radiation force pulse imaging, and divide the input data into a predetermined number matrix
  • a data filtering transformation module is configured to perform forward filtering processing and inverse filtering processing on each data matrix separately, and perform Hilbert transform on the filtered data matrix to obtain analytical data; and a displacement obtaining module, which is used for Performing a cross-correlation operation on the analytical data to obtain a predetermined number of displacement matrices;
  • the displacement combining module is configured to time accumulate each displacement rectangle, take the cumulative and maximum column displacement, and combine the column displacements extracted by the predetermined number of displacement matrices to generate a new displacement matrix;
  • An interpolation module configured to perform cubic spline interpolation in a time direction on the obtained new displacement matrix; a fitting module, configured to perform a Radon transformation on the interpolation result, to obtain a Radon transformation matrix; and an estimation module, configured to acquire the The parameters of the acoustic radiation force pulse imaging are obtained according to the parameters and the Radon transformation matrix.
  • One or more computer storage media containing computer executable instructions for performing an acoustic radiation force pulse imaging estimation method, the method comprising the steps of:
  • the data filtering transformation step performs forward filtering processing and inverse filtering processing on each data matrix, and performs Hilbert transform on the filtered data matrix to obtain analytical data;
  • the displacement combining step performs time accumulation on each displacement rectangle, takes the cumulative and maximum column displacement, and combines the column displacements extracted by the predetermined number of displacement matrices to generate a new displacement matrix;
  • Interpolation step performing cubic spline interpolation in the time direction on the obtained new displacement matrix; fitting step, performing Radon transformation on the interpolation result to obtain a Radon transformation matrix;
  • the estimating step acquires parameters of the acoustic radiation force pulse imaging, and obtains a shear wave rate according to the parameter and the Ladong transform matrix.
  • the above-mentioned acoustic radiation force pulse imaging estimation method and system, computer storage medium by performing forward filtering processing and inverse filtering processing on input data, eliminates relative offset caused by filtering, and performs Hilbert transform on the filtering result, The error generated by the subsequent cross-correlation calculation displacement is made smaller, which makes the calculation result more accurate, and can accurately estimate the small displacement and improve the sensitivity.
  • CUDA parallel computing is used in the estimation process to improve efficiency, and in parallel implementation, it does not increase the amount of repetitive calculation, does not use thread synchronization and shared memory, achieves high thread parallelism, and improves processing efficiency.
  • FIG. 1 is a flow chart of an acoustic radiation force pulse imaging estimation method in an embodiment
  • Figure 2 is a flow chart of the step of determining the displacement
  • Figure 3 is a schematic diagram of pulling transformation
  • FIG. 4 is a structural block diagram of an acoustic radiation force pulse imaging estimation system in one embodiment.
  • Figure 5 is a schematic diagram showing the internal structure of the displacement obtaining module
  • Fig. 6 is a schematic structural view of a terminal for performing an acoustic radiation force pulse imaging estimation method.
  • the acoustic radiation force pulse imaging estimation method can be run on a computer.
  • the sound radiation force pulse imaging estimation method includes:
  • Step 102 The data dividing step reads the input data collected by the acoustic radiation force pulse imaging, and divides the input data into a predetermined number of data matrices.
  • the input data is divided into m data matrices according to the lateral collection position.
  • m is the number of positions in the lateral acquisition position.
  • Each data matrix size is called x col , which means that the row frame data (time direction) is included, and each frame has co/elements (depth direction).
  • Step 104 The data filtering transformation step performs forward filtering processing and inverse filtering processing on each data matrix, and performs Hilbert transform on the filtered data matrix to obtain analytical data.
  • forward filtering processing and inverse filtering processing on the data matrix can eliminate the relative offset caused by the filtering.
  • tOutPut ⁇ k, i) ⁇ (param[j] * input(k, ij)) ( 2 )
  • ram is the filter
  • Z>3 ⁇ 4TO is the matrix to be filtered
  • k respectively represents the filter to be filtered
  • C ⁇ is the filtering result matrix
  • * represents the multiplication.
  • the forward filtering process is processed in parallel using a graphics processing unit (GPU).
  • the input data matrix is the matrix to be filtered.
  • the forward filter parameters are calculated by tools such as Matlab.
  • ⁇ «, >>> is the CUDA syntax symbol, indicating the number of threads that the GPU needs to start the kernel.
  • Input is a structure type, and -> means an element of the structure.
  • the GPU is computationally intensive for data parallel computing.
  • the GPU has more compute cell transistors and processing cores than the CPU (Central Processing Unit), enabling highly parallel computing, such as the Kepler GK110. Consisting of 7.1 billion transistors, it has 15 stream processor cores, each with 196 stream processors, providing over 1 trillion double-precision floating-point calculations per second.
  • CPU Central Processing Unit
  • k respectively represent the first row and the column of the matrix to be filtered.
  • M ⁇ rf is the filtering result matrix of inverse filtering
  • ⁇ is the filter width
  • * means multiplication.
  • the initial content is the same.
  • the inverse filtering process is processed by a graphics processor. Kernel function of the inverse filtering process
  • the filter is the same as the forward filter.
  • input->rows is the number of rows in the input data matrix
  • input->cols is the number of columns in the input data matrix.
  • the input data matrix is the matrix to be filtered.
  • the kernel function Bandpass_back processes an element according to formula (3) and stores the obtained result in the position corresponding to out P ut .
  • Equation (4) the Hilbert transform of the signal r , w is the coordinates of the sample point, t is time, * is a convolution operation.
  • the Fourier transform and the inverse transform are used to implement the Hilbert transform. Specifically, the Fourier transform is performed on the initial data matrix, and the real part of the transform result is multiplied by -1 and then exchanged with the imaginary part. The Fourier normalization inverse transform, and then exchanges the real imaginary part, and then adds the initial matrix to obtain the analytical data.
  • normalization refers to the result of the transformation divided by the number of transformed elements.
  • the Hilbert transform uses a GPU-optimized design, in which both the Fourier transform and the inverse transform use the CUFFT library function cufftExecC2C of CUDA (ie CPU + GPU programming model).
  • the two real imaginary part exchanges and the matrix addition are implemented by two kernel functions.
  • the first kernel function kerne 1_ 1 implements the real imaginary part exchange before the first Fourier transform
  • the second kernel function Kernel_2 implements the second real imaginary part exchange, and adds after the exchange. Since the cufftExecC2C result is not normalized, it is processed during the second real imaginary part exchange.
  • Kernel function call sequence 1) cufftExecC2C (positive transform); 2) kemel l; 3) cufftExecC2C (inverse transform); 4) kernel_2.
  • forward filtering processing, inverse filtering processing, and Hilbert transform processing are sequentially performed on the predetermined number of data matrices in sequence using a plurality of streams.
  • the CPU calls the GPU function interface to transfer the data to be processed to the memory.
  • the GPU reads the data through the memory, calculates the data, and stores the result in the memory.
  • the CPU then calls the data transfer function interface to transfer the memory data back to the memory for the next logical processing.
  • the transfer of data between the CPU and the GPU can be paralleled with the GPU computing process.
  • Data streams can be processed in parallel with the GPU computing process using multiple streams. Stream and stream run in parallel, so although the GPU can only process one kernel function, when multiple streams exist, the GPU can process the kernel function and the memory can transfer data, thus improving the utilization of the device.
  • Step 106 Perform a displacement obtaining step, perform a cross-correlation operation on the parsed data, and obtain a predetermined number of displacement matrices.
  • the displacement obtaining step includes:
  • Step 202 Select a window, and the length of the window is a predetermined value.
  • the window length is M and M is a positive number.
  • Step 204 Calculate the cross-correlation coefficient of each window for the parsed data according to the length of the window, and form a cross-correlation matrix.
  • the window cross-correlation is used to obtain the displacement.
  • the window cross-correlation refers to the cross-correlation only for a segment whose length is the length of the window, and the window position slides with the step length until the end of the data.
  • Equation (5) is the cross-correlation number
  • rA ( w + ", is the analytical data
  • m + u + q, s + i) is the analytical data + ⁇ + ⁇ + complex conjugate
  • s is the sampling time, g For [- ⁇ , ⁇ ].
  • lag ie, the change interval parameter
  • lag takes a range [-N, N], N ⁇ M.
  • WinNum (InputWidth - WindowHW)/Step+l , where WindowHW is the window width, Step is Step, and InputWidth is the number of columns of the input data matrix.
  • WinNum cross-correlation calculation method is as follows: Calculate the relationship between the adjacent rows of Row(i) and Row(i+l) in the data matrix, and then calculate the cross-correlation coefficient of the 0th window.
  • the cross-correlation coefficient of the nth window is calculated as: the cross-correlation coefficient of the n-1th window minus the cross-correlation coefficient of the Step elements before the n-1th window, plus the mutual interaction of the last Step elements of the nth window The number of relationships.
  • the number of cross-correlations of the first window is: the number of cross-correlations of the 0th window minus the number of cross-correlations of the Step elements before the 0th window, plus the cross-correlation of the last Step elements of the first window;
  • the number of cross-correlations of the second window the number of cross-correlations of the first window minus the number of cross-correlations of the Step elements in the first window, plus the number of cross-correlations of the last Step elements of the second window, calculated in turn Until the last window.
  • the kernel function ComputeXcross calculates a lag correlation coefficient.
  • Step 206 Obtain the change interval parameter value corresponding to the maximum cross-correlation coefficient in the cross-correlation matrix.
  • the kernel function MatDelay is used to find the extremum and displacement of the cross-correlation matrix. Find the lag value lag_ corresponding to the largest correlation coefficient in ⁇ ′′
  • Step 208 Determine a displacement according to the change interval parameter value.
  • Delay ( lag _ max + delta _ t ) * factor ( 6 ) where delta_t is obtained using equation (7).
  • Fa ct. r c / ( 2 *fe), c is the sound speed 1540m / s , fs sampling frequency, * means multiplication.
  • the mutual relation number can be obtained by using the CUDA flow method.
  • the CUDA stream can hide data transfer time and further reduce program run time.
  • the stream processing example for the displacement is as follows:
  • ⁇ Data is input into memory by memory.
  • the forward filtering process, the inverse filtering process, and the Hilbert transform can be performed by means of the stream.
  • Step 108 The displacement combining step performs time accumulation on each displacement rectangle, takes the cumulative and maximum column displacement, and combines the column displacements extracted by the predetermined number of displacement matrices to generate a new displacement matrix.
  • the obtained m displacement matrices are respectively processed: performing time (column direction) accumulation, taking cumulative and maximum column displacements, and combining the column displacements extracted by the m displacement matrices to generate a new displacement matrix, ie, m X (row - 1).
  • Step 110 The interpolation step performs cubic spline interpolation in the time direction on the obtained new displacement matrix.
  • the interpolation step comprises: using a central processing unit to obtain coefficients of a cubic spline interpolation function between each interpolation segment, and then using a graphics processor to calculate an interpolation multiple of each segment and subtracting one interpolation point according to the cubic spline interpolation function. Value, where the number of blocks of the kernel function called by the graphics processor The quantity is the number of rows of the matrix to be interpolated, and the number of threads is the number of columns of the matrix to be interpolated.
  • the cubic spline interpolation multiple can be set as needed.
  • the linear equation matrix of the interpolation coefficient is a three-diagonal matrix, which is solved by direct LU decomposition. Since the direct decomposition of LU is an iterative process, the CPU calculates all the coefficients of the interpolation function, and then uses the GPU to calculate the values of all the interpolation points in parallel.
  • spline_kern Design the kernel function spline_kern to calculate all the values of a [ ⁇ ⁇ " ] interpolation segment.
  • Step 112 a fitting step, performing a Radon transformation on the interpolation result to obtain a Radon transformation matrix.
  • the pull map is pulled.
  • the sum of the points from the point to the opposite point is sequentially determined.
  • the sequence number of the end point cannot be smaller than the sequence number of the starting point.
  • the matrix width is Q
  • the dimension of the grid and the block are usually directly related to the length and width of the matrix.
  • the value of Q is usually larger than 1960, then the size of the Radon transformation matrix is 1960*1960.
  • the maximum number of threads in the GPU supporting block with a computing power of 3.5 is 1024, so the number of internal threads in the block cannot be defined by Q.
  • the x-dimension of the grid and the x-dimension of the block represent the serial number of the starting endpoint index_start (0 to 1959), and the y-dimension of the grid and the y-dimension of the block represent the serial number index_end of the end point.
  • the kernel function implements a diagonal diagonal summation (summing the data from the initial point index_start to the diagonal of the end point index_end), and the resulting result Sum (index start, index end ) is saved to the output data.
  • Step 114 The estimating step acquires a parameter of the sound radiation force pulse imaging, and obtains a shear wave rate according to the parameter and the Radon transformation matrix.
  • the parameters of the acoustic radiation force pulse imaging in the estimating step include the spacing of the two array elements in the probe, the number of probes, the number of rows of the input data matrix, the frame rate, and the sampling multiple.
  • Vs pitch * pitchNum * (dataRow - ⁇ ) I ⁇ L * fs * spline)
  • the above-mentioned acoustic radiation force pulse imaging estimation method eliminates the relative offset caused by the filtering by performing forward filtering processing and inverse filtering processing on the input data, and performs Hilbert transform on the filtering result, thereby reducing subsequent cross-correlation calculation.
  • the error generated during the displacement makes the calculation result more accurate, and the micro displacement can be accurately estimated, and the sensitivity is improved.
  • FIG. 4 it is a structural block diagram of an acoustic radiation force pulse imaging estimation system in one embodiment. Details of the acoustic radiation force pulse imaging estimation system that are not described in detail can be referred to the corresponding method embodiments described above.
  • the acoustic radiation force pulse imaging estimation system comprises a data dividing module 410, a data filtering transformation module 420, a displacement obtaining module 430, a displacement combining module 440, an interpolation module 450, a fitting module 460, and an estimating module 470. among them:
  • the data dividing module 410 is configured to read the input data acquired by the acoustic radiation force pulse imaging, and divide the input data into a predetermined number matrix. Specifically, the input data is divided into m data matrices according to the lateral gather position. m is the number of positions in the lateral collection position. Each data matrix size is r OW X CO /, which means that it contains row frame data (time direction), and each frame has co/elements (depth direction).
  • the data filtering transformation module 420 is configured to perform forward filtering processing and inverse filtering processing on each data matrix respectively, and perform Hilbert transform on the filtered data matrix to obtain analytical data.
  • tOutPutik, ⁇ ) ⁇ (param[j] * input k, i - j)) where param m is the filter length, /3 ⁇ 4rot is the matrix to be filtered, and k is the row and column of the matrix to be filtered, respectively.
  • OWPW is the filtering result matrix, and * means multiplication.
  • the input data matrix is the matrix to be filtered.
  • k respectively represent the first row and the column of the matrix to be filtered.
  • i ⁇ is the filtered result matrix of the inverse filtering, which is the filter width, and * represents the multiplication.
  • the inverse filtering process is processed by a graphics processor. Kernel function of the inverse filtering process
  • Bandpass back (the kernel function called by the graphics processor)
  • the kernel function Bandpass_back processes an element according to formula (3) and stores the obtained result in the position corresponding to out P ut .
  • the data filter transform module 420 performs the Hilbert transform formula as:
  • ⁇ ( " ⁇ ) is a signal r (w Hilbert transform, w is The coordinates of the sampling point are time, and * is a convolution operation.
  • the Fourier transform and the inverse transform are used to implement the Hilbert transform. Specifically, the Fourier transform is performed on the initial matrix, and the real part of the transform result is multiplied by -1 and then exchanged with the imaginary part. The inverse transformation of the inner leaf is performed, and the real imaginary part is exchanged, and then the initial matrix is added to obtain the analytical data.
  • the Hilbert transform uses a GPU-optimized design, in which both the Fourier transform and the inverse transform use the CUFFT library function cufftExecC2C of CUDA (ie CPU + GPU programming model).
  • the two real imaginary part exchanges and the matrix addition are implemented by two kernel functions.
  • the first kernel function kerne 1_ 1 implements the real imaginary part exchange before the first Fourier transform
  • the second kernel function Kernel_2 implements the second real imaginary part exchange, and adds after the exchange. Since the cufftExecC2C result is not normalized, it is processed during the second real imaginary part exchange.
  • Kernel function call sequence 1) cufftExecC2C (positive transform); 2) kemel l; 3) cufftExecC2C (inverse transform); 4) kernel_2.
  • the data filtering transformation module 420 further performs forward filtering processing, inverse filtering processing, and Hilbert transform processing on the predetermined number of data matrices in sequence using a plurality of streams.
  • one of the characteristics of GPU computing is that the data to be processed needs to be transferred from the memory to the memory, and the result data needs to be transferred from the memory to the memory.
  • the CPU calls the GPU function interface to transfer the data to be processed to the memory.
  • the GPU reads the data through the memory, calculates the data, and stores the result.
  • the memory is stored, and the CPU calls the data transfer function interface to transfer the memory data back to the memory for the next logical processing.
  • the transfer of data between the CPU and the GPU can be parallel to the GPU computing process.
  • Data streams can be processed in parallel with the GPU computing process using multiple streams.
  • the flow and the flow run in parallel, so although the GPU can only process one kernel function, when multiple streams exist, the GPU can process the kernel function and the memory can transmit data, thereby improving the utilization of the device.
  • the displacement obtaining module 430 is configured to perform a cross-correlation operation on the analytical data to obtain a predetermined number of displacement matrices.
  • the displacement obtaining module 430 includes a selection sub-module 432, a cross-correlation calculation sub-module 434, a comparison sub-module 436, and a displacement calculation sub-module 438. among them:
  • the selection sub-module 432 is used to select a window whose length is a predetermined value.
  • the window length is M and M is a positive number.
  • the cross-correlation calculation sub-module 434 is configured to calculate the cross-correlation coefficient of each window for the parsed data according to the length of the window, and form a cross-correlation matrix.
  • the window cross-correlation is used to obtain the displacement.
  • the window cross-correlation refers to the cross-correlation only for a segment whose length is the length of the window, and the window position slides with the step length until the end of the data.
  • Equation (5) is the cross-correlation coefficient, + M , is the analytical data , rm + u + q, s + ⁇ ) is the analytical data + ⁇ + ⁇ + complex conjugate, s is the sampling time, g is [- ⁇ , ⁇ ].
  • lag ie, the change interval parameter
  • lag takes a range [-N, N], N ⁇ M.
  • WinNum (InputWidth - WindowHW)/Step+l , where WindowHW is the window width, Step is Step, and InputWidth is the number of columns of the input data matrix.
  • WinNum cross-correlation calculation method is as follows: Calculate the relationship between the adjacent rows of Row(i) and Row(i+l) in the data matrix, and then calculate the cross-correlation coefficient of the 0th window.
  • the cross-correlation coefficient of the nth window is calculated as: the cross-correlation coefficient of the n-1th window minus the cross-correlation coefficient of the Step elements before the n-1th window, plus the mutual interaction of the last Step elements of the nth window Relationship Number.
  • the number of cross-correlations of the first window is: the number of cross-correlations of the 0th window minus the number of cross-correlations of the Step elements before the 0th window, plus the cross-correlation of the last Step elements of the first window;
  • the number of cross-correlations of the second window the number of cross-correlations of the first window minus the number of cross-correlations of the Step elements in the first window, plus the number of cross-correlations of the last Step elements of the second window, calculated in turn Until the last window.
  • the kernel function ComputeXcross calculates a lag mutual relation number, and the number of threads designing one block is equal to the total number of lags, and the correlation coefficient of all windows in the adjacent two rows can be calculated.
  • the matrix needs to be calculated (rows-1) and the windowed cross-correlation (rows is the number of rows in the matrix), which is implemented by (rows -1) blocks.
  • a window's cross-correlation coefficient (with lagNum cross-correlation coefficient) is represented by ⁇ ′′.
  • the WinNum-1 comparison sub-module 436 is configured to obtain the change interval parameter value corresponding to the maximum cross-correlation coefficient in the cross-correlation matrix.
  • the kernel function MatDelay is used to find the extremum and displacement of the cross-correlation matrix. Find the lag value lag_max corresponding to the largest correlation coefficient in ⁇ ′′
  • the displacement calculation sub-module 438 is configured to determine the displacement based on the change interval parameter value.
  • Delay ( lag _ max + delta _ t ) * factor ( 6 ) where delta_t is obtained using equation (7).
  • ct. r c / ( 2 *fe), c is the sound speed 1540m / s , fs sampling frequency, * means multiplication.
  • the small displacement d(m, S ) can be obtained by the ultrasonic speed C :
  • the displacement combining module 440 is configured to time accumulate each displacement rectangle, take the cumulative and maximum column displacement, and combine the column displacements extracted by the predetermined number of displacement matrices to generate a new displacement matrix. Specifically, the obtained m displacement matrices are respectively processed: performing time (column direction) accumulation, taking cumulative and maximum column displacements, and combining the column displacements extracted by the m displacement matrices to generate a new displacement matrix, ie, mx 1).
  • the interpolation module 450 is configured to perform cubic spline interpolation in the time direction on the obtained new displacement matrix.
  • the interpolation module 450 is further configured to use a central processor to obtain coefficients of a cubic spline interpolation function between each interpolation segment, and then use a graphics processor to calculate each segment in parallel according to the cubic spline interpolation function.
  • the interpolation multiple is reduced by the value of one interpolation point, wherein the number of blocks of the kernel function called by the graphics processor is the number of rows of the matrix to be interpolated, and the number of threads is the number of columns of the matrix to be interpolated.
  • a ', b are interpolation function coefficients, which can be solved by linear Equation is obtained, Since the linear equation is an indefinite solution equation, the natural solution can be used to determine the equation solution.
  • the linear equation matrix of the interpolation coefficient is a three-diagonal matrix, which is solved by direct LU decomposition. Since the direct decomposition of LU is an iterative process, the CPU calculates all the coefficients of the interpolation function, and then uses the GPU to calculate the values of all the interpolation points in parallel.
  • spline_kern Design the kernel function spline_kern to calculate all the values of a [ ⁇ ⁇ " ] interpolation segment.
  • the fitting module 460 is configured to perform a Radon transform on the interpolation result to obtain a Radon transformation matrix.
  • the pull map is pulled.
  • the sum of the points from the point to the opposite point is sequentially determined.
  • the sequence number of the end point cannot be smaller than the sequence number of the starting point.
  • the matrix width is Q
  • the dimension of the grid and the block are usually directly related to the length and width of the matrix.
  • the value of Q is usually larger than 1960, then the size of the Radon transformation matrix is 1960*1960.
  • the maximum number of threads in the GPU supporting block with a computing power of 3.5 is 1024, so the number of internal threads in the block cannot be defined by Q.
  • the x-dimension of the grid and the x-dimension combination of the block represent the index of the starting endpoint index start (0 to 1959), and the y-dimension of the grid and the y-dimension of the block represent the index of the end point index_end.
  • the kernel function implements a diagonal diagonal summation (summing the data from the initial point index_start to the diagonal of the end point index_end), and the result Sum ( index start, index end ) is saved to the output data.
  • the estimating module 470 is configured to obtain a parameter of the acoustic radiation force pulse imaging, and obtain a shear wave rate according to the parameter and the Radon transformation matrix.
  • the parameters of the acoustic radiation force pulse imaging include the spacing of two array elements in the probe, The number of probes, the number of rows of the input data matrix, the frame rate, and the multiple of the sample;
  • the formula for the estimation module to determine the shear wave rate is:
  • Vs pitch * pitchNum * (dataRow - ⁇ ) I ⁇ L * fs * spline)
  • the above-mentioned acoustic radiation force pulse imaging estimation system eliminates the relative offset caused by the filtering by performing forward filtering processing and inverse filtering processing on the input data, and performs Hilbert transform on the filtering result, thereby reducing subsequent cross-correlation calculation.
  • the error generated during the displacement makes the calculation result more accurate, and the micro displacement can be accurately estimated, and the sensitivity is improved.
  • CUDA parallel computing is used in the estimation process to improve efficiency, and in parallel implementation, it does not increase the amount of repetitive calculation, does not use thread synchronization and shared memory, achieves high thread parallelism, and improves processing efficiency.
  • This estimation method uses 4 sets of 100*512 data as a set of input data, and sets the window length 100, the total number of lags is 81, and the time of calculation is 80ms.
  • This time-consuming time includes the data transfer time between the GPU and the CPU and the GPU calculation time, excluding the GPU application memory time. Since the GPU application for video memory is slow (approx. 500ms for applying enough memory), considering that the applied memory can be reused, you can apply for enough memory at the system startup time and release it when the system is shut down. The application memory time is not calculated as the running time. It is reasonable. Before the GPU optimization, the algorithm running on a personal computer configured with dual-core 3.3G and 3.29G memory 2G takes about 2000ms.
  • FIG. 6 is a schematic structural diagram of a terminal for operating the acoustic radiation force pulse imaging estimation method and system according to the present invention.
  • the terminal can be used to implement the acoustic radiation force pulse imaging estimation method provided in the above embodiment. Specifically:
  • Terminal 1200 can include a memory having one or more computer readable storage media 110.
  • the terminal structure shown in FIG. 6 does not constitute a limitation to the terminal, and may include more or less components than those illustrated, or a combination of certain components, or different component arrangements. among them:
  • the memory 110 can be used to store software programs and modules, and the processor 150 executes various functional applications and data processing by running software programs and modules stored in the memory 110.
  • the memory 110 may mainly include a storage program area and a storage data area, wherein the storage program area may store an operating system, an application required for at least one function (such as a sound playing function, an image playing function, etc.), and the like; the storage data area may be stored according to The data created by the use of the terminal 1200 (such as acoustic radiation force pulse imaging data) or the like.
  • memory 110 can include high speed random access memory, and can also include non-volatile memory, such as at least one magnetic disk storage device, flash memory device, or other volatile solid state storage device. Accordingly, memory 110 may also include a memory controller to provide access to memory 110 by processor 150 and input unit 120.
  • Input unit 120 can be used to receive input numeric or character information, as well as to generate keyboard, mouse, joystick, optical or trackball signal inputs related to user settings and function controls.
  • input unit 120 can include touch-sensitive surface 121 as well as other input devices 122.
  • the touch-sensitive surface 121 also referred to as a touch display or trackpad, can collect touch operations on or near the user (eg, the user uses a finger, stylus, etc., on any suitable object or accessory on the touch-sensitive surface 121 or The operation near the touch-sensitive surface 121) and driving the corresponding connecting device according to a preset program.
  • the touch-sensitive surface 121 can include two portions of a touch detection device and a touch controller.
  • the touch detection device detects the touch orientation of the user, and detects a signal brought by the touch operation, and transmits the signal to the touch controller; the touch controller receives the touch information from the touch detection device, converts the touch information into contact coordinates, and sends the touch information
  • the processor 150 is provided and can receive commands from the processor 150 and execute them.
  • the touch-sensitive surface 121 can be realized by various types such as resistive, capacitive, infrared, and surface acoustic waves.
  • the input unit 120 can also include other input devices 122.
  • other input devices 122 may include, but are not limited to, one or more of a physical keyboard, function keys (such as volume control buttons, switch buttons, etc.), trackballs, mice, joysticks, and the like.
  • the display unit 130 can be used to display information input by the user or information provided to the user and various graphical user interfaces of the terminal 1200, which can be composed of graphics, text, icons, video, and any combination thereof.
  • the display unit 130 may include a display panel 131, optionally,
  • the display panel 131 is configured in the form of an LCD (Liquid Crystal Display), an OLED (Organic Light-Emitting Diode), or the like.
  • the touch-sensitive surface 121 can cover the display panel 131, and when the touch-sensitive surface 121 detects a touch operation thereon or nearby, it is transmitted to the processor 150 to determine the type of the touch event, and then the processor 150 according to the touch event The type provides a corresponding visual output on display panel 131.
  • touch-sensitive surface 121 and display panel 131 are implemented as two separate components to implement input and input functions, in some embodiments, touch-sensitive surface 121 can be integrated with display panel 131 for input. And output function.
  • the audio circuit 140, the speaker 141, and the microphone 142 provide an audio interface between the user and the terminal 1200.
  • the audio circuit 140 can transmit the converted electrical data of the received audio data to the speaker 141 for conversion to the sound signal output by the speaker 141; on the other hand, the microphone 142 converts the collected sound signal into an electrical signal, and the audio circuit 140 After receiving, it is converted into audio data, and then processed by the audio data output processor 150.
  • the audio circuit 140 may also include an earbud jack to provide communication of the peripheral earphones with the terminal 1200.
  • the processor 150 is a control center of the terminal 1200 that connects various portions of the entire computer using various interfaces and lines, by running or executing software programs and/or modules stored in the memory 110, and recalling data stored in the memory 110, Various functions and processing data of the terminal 1200 are executed.
  • the processor 150 may include one or more processing cores.
  • the processor 150 may integrate an application processor and a modem processor, where the application processor mainly processes an operating system, a user interface, an application, and the like.
  • the modem processor primarily handles wireless communications. It will be appreciated that the above described modem processor may also not be integrated into the processor 150.
  • the terminal 1200 also includes a power source 160 (such as a battery) for powering various components.
  • the power source can be logically coupled to the processor 150 through a power management system to implement functions such as managing charging, discharging, and power management through the power management system.
  • the power supply 160 can also include any one or more of a DC or AC power source, a recharging system, a power failure detection circuit, a power converter or an inverter, a power status indicator, and the like.
  • the terminal 1200 may further include a camera, a Bluetooth module, and the like, and details are not described herein.
  • the terminal 1200 is a computer.
  • the invention also provides one or more computer storage media comprising computer executable instructions comprising one or more programs, the one or more programs being executed by a processor of one or more terminals
  • An acoustic radiation force pulse imaging estimation method includes the following steps:
  • the data filtering transformation step performs forward filtering processing and inverse filtering processing on each data matrix, and performs Hilbert transform on the filtered data matrix to obtain analytical data;
  • the displacement combining step performs time accumulation on each displacement rectangle, takes the cumulative and maximum column displacement, and combines the column displacements extracted by the predetermined number of displacement matrices to generate a new displacement matrix;
  • Interpolation step performing cubic spline interpolation in the time direction on the obtained new displacement matrix; fitting step, performing Radon transformation on the interpolation result to obtain a Radon transformation matrix;
  • the estimating step acquires parameters of the acoustic radiation force pulse imaging, and obtains a shear wave rate according to the parameter and the Ladong transform matrix.
  • forward filtering processing, inverse filtering processing, and Hilbert transform processing are performed on the predetermined number of data matrices in sequence by using a plurality of streams.
  • the formula of the forward filtering process in the data filtering transformation step is: when i ⁇ iParalen:
  • Output ( i, Width - 1 - /) ⁇ (param[j] * input ⁇ k, Width _ 1 + j— /))
  • k respectively represent the first row and the column of the matrix to be filtered.
  • i ⁇ is the filtered result matrix of the inverse filtering, which is the filter width, and * represents the multiplication;
  • the kernel function and the inverse filtering processing of the forward filtering process are both processed in parallel by a graphics processor, and the number of blocks of the kernel function called by the graphics processor is the same as the number of threads, and the number of the blocks is the number of rows of the matrix to be filtered.
  • the number of threads is the number of columns of the matrix to be filtered.
  • the Hilbert transform in the data filtering transformation step comprises: performing a Fourier transform on the initial data matrix, multiplying the real part of the transformation result by -1, and exchanging with the imaginary part, and then performing the Fu The inverse transformation is performed on the inner leaf, and the real imaginary part is exchanged, and then the initial data matrix is added to obtain the analytical data.
  • the step of determining the displacement comprises:
  • the interpolating step comprises:
  • the central processor is used to obtain the coefficients of the cubic spline interpolation function between each interpolation segment, and then the graphics processor is used to calculate the interpolation multiple of each segment and the value of one interpolation point in parallel according to the cubic spline interpolation function, wherein
  • the number of blocks of the kernel function called by the graphics processor is the number of rows of the matrix to be interpolated, and the number of threads is the number of columns of the matrix to be interpolated.
  • the parameters of the acoustic radiation force pulse imaging in the estimating step include a spacing of two array elements in the probe, a number of probes, a number of rows of the input data matrix, a frame rate, and a sampling multiple;
  • the formula for obtaining the shear wave velocity in the estimating step is:
  • Vs pitch * pitchNum * (dataRow - ⁇ ) I ⁇ L * fs * spline)
  • the storage medium may be a magnetic disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM). It is not to be construed as limiting the scope of the invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the spirit and scope of the invention. Therefore, the scope of the invention is to be determined by the appended claims.

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Pathology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Biophysics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Remote Sensing (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种声辐射力脉冲成像估算方法和系统、计算机存储介质。所述方法包括：读取输入数据，并将所述输入数据分成预定数量数据矩阵；分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；对所述解析数据进行互相关运算，求取得到预定数量位移矩阵；对每个位移矩形累加，取累积和最大的列位移，生成新的位移矩阵；对得到的新的位移矩阵进行时间方向上的三次样条插值；对插值结果进行拉东变换，得到拉东变换矩阵；获取所述声辐射力脉冲成像的参数，根据所述参数及拉东变换矩阵求取剪切波速率。减小了互相关计算位移时产生的误差，计算结果更加准确，且对微小位移能精确估算，提高了灵敏度。

Description

说 明 书
发明名称： 声辐射力脉冲成像估算方法和系统、 计算机存储介盾
【技术领域】
本发明涉及医学影像领域， 特别是涉及一种声辐射力脉冲成像估算方法 和系统、 计算机存储介质。
【背景技术】
声辐射力脉冲成像 ( Acoustic Radiation Force Impulse, ARFI )是通过利 用医学超声功率范围内的聚焦超声波束产生的辐射力， 使生物粘弹性组织局 部区域产生微小变形并形成沿侧向传播的剪切波， 利用超声成像技术对微小 形变进行追踪并计算出剪切波的侧向传播速度， 最后利用弹性重构算法由剪 切波速度估算出组织的弹性分布的一种成像技术。
声辐射力脉冲成像计算是指估算组织内部剪切波的位移， 后期通过对位 移的判断了解组织情况。位移估算方法有多种，例如基于通用 GPU ( Graphics Process Unit, 图形处理器）编程的实时超声成像监测方法， 步骤包括： 1 )展 开数据； 2 )加窗数据； 3 ) 互相关处理数据； 4 )估计位移； 5 )在整个时间 深度上重复步骤 3和 4; 6 ) 中值滤波； 7 )估计速率。 该监测方法使用了普 通的基于时域互相关算法， 针对微小位移的估算灵敏度不高， 准确度低。
【发明内容】
基于此， 有必要针对微小位移估算灵敏度低的问题， 提供一种灵敏度高 的声辐射力脉冲成像估算方法。
此外， 还有必要提供一种灵敏度高的声辐射力脉冲成像估算系统。
一种声辐射力脉冲成像估算方法， 包括：
数据划分步骤， 读取通过声辐射力脉冲成像所采集的输入数据， 并将所 述输入数据分成预定数量数据矩阵； 数据滤波变换步骤， 分别对每个数据矩阵进行前向滤波处理和逆向滤波 处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；
位移求取步骤， 对所述解析数据进行互相关运算， 求取得到预定数量位 移矩阵；
位移合并步骤， 对每个位移矩形进行时间累加， 取累积和最大的列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值步骤， 对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合步骤， 对插值结果进行拉东变换， 得到拉东变换矩阵；
估算步骤， 获取所述声辐射力脉冲成像的参数， 根据所述参数及拉东变 换矩阵求取剪切波速率。
一种声辐射力脉冲成像估算系统， 包括：
数据划分模块， 用于读取通过声辐射力脉冲成像所釆集的输入数据， 并 将所述输入数据分成预定数量矩阵；
数据滤波变换模块， 用于分别对每个数据矩阵进行前向滤波处理和逆向 滤波处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据； 位移求取模块， 用于对所述解析数据进行互相关运算， 求取得到预定数 量位移矩阵；
位移合并模块， 用于对每个位移矩形进行时间累加， 取累积和最大的列 位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值模块，用于对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合模块， 用于对插值结果进行拉东变换， 得到拉东变换矩阵； 估算模块， 用于获取所述声辐射力脉冲成像的参数， 根据所述参数及拉 东变换矩阵求取剪切波速率。
此外， 还提供一种计算机存储介质。
一个或多个包含计算机可执行指令的计算机存储介质， 所述计算机可执 行指令用于执行一种声辐射力脉冲成像估算方法， 其特征在于， 所述方法包 括以下步骤：
数据划分步骤， 读取通过声辐射力脉冲成像所釆集的输入数据， 并将所 述输入数据分成预定数量数据矩阵；
数据滤波变换步骤， 分别对每个数据矩阵进行前向滤波处理和逆向滤波 处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；
位移求取步骤， 对所述解析数据进行互相关运算， 求取得到预定数量位 移矩阵；
位移合并步骤， 对每个位移矩形进行时间累加， 取累积和最大的列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值步骤， 对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合步骤， 对插值结果进行拉东变换， 得到拉东变换矩阵；
估算步骤， 获取所述声辐射力脉冲成像的参数， 根据所述参数及拉东变 换矩阵求取剪切波速率。
上述声辐射力脉冲成像估算方法和系统、 计算机存储介质， 通过对输入 数据进行前向滤波处理和逆向滤波处理， 消除了滤波造成的相对偏移， 并对 滤波结果进行希尔伯特变换， 减小了后续互相关计算位移时产生的误差， 使 得计算结果更加准确， 且对微小位移能精确估算， 提高了灵敏度。
另外， 估算过程中釆用 CUDA并行计算实现， 提高了效率， 且在并行实 现上， 不增加重复计算量， 不使用线程同步和共享内存， 实现线程高度并行， 提高了处理效率。
【附图说明】
图 1为一个实施例中声辐射力脉冲成像估算方法的流程图；
图 2为该位移求取步骤的流程图；
图 3为拉动变换示意图；
图 4为一个实施例中声辐射力脉冲成像估算系统的结构框图。
图 5为该位移求取模块的内部结构示意图；
图 6为执行声辐射力脉冲成像估算方法的终端的结构示意图。
【具体实施方式】 为了使本发明的目的、 技术方案及优点更加清楚明白， 以下结合附图及 实施例， 对本发明进行进一步详细说明。 应当理解， 此处所描述的具体实施 例仅仅用以解释本发明， 并不用于限定本发明。
如图 1所示， 为一个实施例中声辐射力脉冲成像估算方法的流程图。 该 声辐射力脉冲成像估算方法可运行于计算机上。 该声辐射力脉冲成像估算方 法， 包括：
步骤 102, 数据划分步骤， 读取通过声辐射力脉冲成像所采集的输入数 据， 并将所述输入数据分成预定数量数据矩阵。
具体的， 将输入数据按侧向采集位置分为 m个数据矩阵。 m为侧向采集 位置的条数。每个数据矩阵大小为謂 x col ,表示包含 row帧数据 (时间方向）， 每帧有 co/个元素 （深度方向）。
步骤 104, 数据滤波变换步骤， 分别对每个数据矩阵进行前向滤波处理 和逆向滤波处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数 据。
具体的， 对数据矩阵进行前向滤波处理和逆向滤波处理， 可消除滤波造 成的相对偏移。
本实施例中， 1 )该数据滤波变换步骤中前向滤波处理的公式为： 当 i < iParalen时：
tOutPut{k, ΐ) = ^ (param[j] * input(k, i j)) + ^ (param[j] * input kfi)) ( 1 ) 当 〉 iParalen时：
tOutPut{k, i) = ^ (param[j] * input(k, i j)) ( 2 ) 其中， ram是滤波器, 是滤波器长度， Z>¾TO为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， C ^是滤波结果矩阵， *表示乘法。
该前向滤波处理采用图形处理器（GPU, Graphics Process Unit ) 并行处 理。 该前向滤波处理的内核函数 Bandpass_front (即图形处理器调用的内核函 数)调用的线程数设计： Bandpass_front«< blocks, threads»>( input, aram, iParaLen , tOutPut) , 此处 blocks = dim3 (input->rows) ； threads = dim3(input->cols); 其中， input ->rows是输入数据矩阵的行数， input->cols是 输入数据矩阵的列数。 该输入数据矩阵即为待滤波矩阵。 前向滤波器参数由 Matlab等工具计算求得。 其中， <«, >>>是 CUDA语法符号， 表示 GPU执 行内核需要启动的线程数， 使用方式： <<<表示启动的 block数量， 表示 1个 block里面的线程数 >>>。 input为一个结构类型， ->表示取结构的一个元素。
GPU为面向计算密集型处理，适用数据并行化计算， GPU与 CPU( Central Processing Unit, 中央处理器）相比具有更多的计算单元晶体管和处理核心， 可实现高度并行化计算， 如 Kepler GK110由 71亿个晶体管组成， 拥有 15个 流处理器核心， 每个核心有 196个流处理器， 提供每秒超过 1万亿次双精度 浮点计算的吞吐量。
2 )该数据滤波变换步骤中逆向滤波处理的公式为：
当 i > iParalen时： output(k, Width - \ - ΐ) = ^ (param[j] * input(k, Width - 1 + j - i ) ( 3 )
7=0
其中， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， 。M^rf为逆向滤波的滤波结果矩阵， ^^为滤波器宽度， *表示乘法。 。 的 初始内容与 相同。
该逆向滤波处理采用图形处理器处理。 该逆向滤波处理的内核函数
Bandpass back (即图形处理器调用的内核函数） 调用的线程数设计： Bandpass_back«< blocks , threads»>() , 其中, blocks = dim3 (input ->rows); threads = dim3 (input->cols) , 滤波器和前向滤波器相同， input->rows是输入数 据矩阵的行数， input->cols是输入数据矩阵的列数。 该输入数据矩阵即为待滤 波矩阵。 内核函数 Bandpass_back将一个元素按照公式 （3 ) 处理， 将获取的 结果存放在 outPut对应的位置。
3 )该数据滤波变换步骤中希尔伯特变换的公式为：
rA (m, t) = r(m, t) * (\ Ι(τα)) ( 4 ) 式（4 ) 中， 为信号 r 的希尔伯特变换， w为采样点的坐标， t是时间， *是卷积运算。
在算法上采用傅里叶变换和逆变换来实现希尔伯特变换， 具体为： 对初 始的数据矩阵进行傅里叶变换， 将变换结果的实部乘 -1后和虚部交换， 再进 行傅里叶归一化逆变换， 再交换实虚部， 再和初始的矩阵相加得到解析数据。 其中， 归一化是指变换结果除以变换元素个数。
希尔伯特变换采用 GPU优化设计， 其中， 傅里叶变换和逆变换均使用 CUDA (即 CPU+GPU编程模型 ) 的 CUFFT库函数 cufftExecC2C。
其中， 两次实虚部交换以及矩阵相加用 2个内核函数来实现， 第一个内核 函数 kerne 1_ 1实现第一次的傅里叶变换前的实部虚部交换， 第二个内核函数 kernel_2实现第二次实部虚部交换， 并且交换后相加， 因 cufftExecC2C结果没 有归一化处理， 在第二次实虚部交换时候处理。
内核函数调用顺序： 1 ) cufftExecC2C (正变换）； 2 ) kemel l； 3 ) cufftExecC2C (逆变换）； 4 ) kernel_2。
两次内核函数的 blocks数量和 threads数量相同，均为 blocks = dim3 (input -〉 rows); threads = dim3 (input -〉 cols), input -〉 rows是输入数据矩阵的行数, input ->cols是输入数据矩阵的列数。
在一个实施例中， 该数据滤波变换步骤中， 采用多个流按顺序依次对预 定数量数据矩阵进行前向滤波处理、 逆向滤波处理和希尔伯特变换处理。 具 体的， 因 GPU计算的一个特点是待处理数据需要从内存传送到显存， 结果数 据需要从显存传回内存这两个数据传输过程。 由 CPU调用 GPU函数接口把 待处理数据传输到显存， GPU通过显存读取数据， 计算数据， 把结果存放显 存， CPU再调用数据传输函数接口将显存数据传回内存， 进行下一步逻辑处 理。 数据在 CPU和 GPU之间的传输过程可与 GPU计算过程并行。 可使用多 个流（stream ) 实现数据传输与 GPU计算过程并行处理。 流与流之间是并行 运行的， 因此虽然 GPU —次只能处理一个内核函数， 而当多个流存在时， GPU在处理内核函数的同时， 显存可以传输数据， 从而提高设备的利用率。
步骤 106, 位移求取步骤， 对该解析数据进行互相关运算， 求取得到预 定数量位移矩阵。 如图 2所示， 该位移求取步骤包括：
步骤 202, 选取窗口， 该窗口长度为预定值。
具体的， 需要选取合适的窗口， 窗口太小会导致计算距离不准确， 太大 会增加不必要的计算量。 本实施例中， 窗口长度为 M , M为正数。
步骤 204, 根据该窗口长度按步长滑动的方式对该解析数据计算每个窗 口的互相关系数， 并形成互相关系数矩阵。
本实施例中， 采用加窗互相关求取位移， 加窗互相关是指一次只针对长 度为窗口长度的一段求互相关， 窗口位置随步长滑动， 直到数据末位为止。
每个窗口的互相关系数计算公式为：
M 12-1
k q, s) = ^ rA (m + u, s)rA (m + u + q, s + \) ( 5 )
u=M / 2
式（5 )中， 是互相关系数， rA(w + "， 是解析数据， m + u + q, s + i) 是解析数据 + Μ + ^ + 求复共轭， s是采样时间， g为 [-Μ,Μ]。
g每一次变化称一个 lag (即变化间隔参数）， lag的范围 [-Μ,Μ]。 当估算 位移比较小时， lag取一个范围 [-N, N] , N < M。
采用内核函数 ComputeXcross计算所有窗（ WinNum个 )中一次 lag的互 相关系数。 WinNum=(InputWidth - WindowHW)/Step+l , 其中， WindowHW 是窗口宽度， 步长为 Step, InputWidth为输入数据矩阵的列数， 为了减少重 复计算提高时间效率， WinNum个互相关系数计算方法如下： 计算数据矩阵 相邻两行 Row(i), Row(i+l)互相关系数， 则先计算第 0个窗口的互相关系数。 第 n个窗口的互相关系数计算为：第 n-1个窗口的互相关系数减去第 n-1个窗 口前 Step个元素的互相关系数，加上第 n个窗口最后 Step个元素的互相关系 数。 例如， 第一个窗口的互相关系数为： 第 0个窗口的互相关系数减去第 0 个窗口前 Step个元素的互相关系数，加上第一个窗口最后 Step个元素的互相 关系数； 第二个窗口的互相关系数： 第一个窗口的互相关系数减去第一个窗 前 Step个元素的互相关系数，加上第二个窗口最后 Step个元素的互相关系数， 依次计算到最后一个窗口为止。
内核函数 ComputeXcross计算 1个 lag互相关系数， 设计 1个 block的线程数 等于 lag总数， 可以计算出相邻两行所有窗口的互相关系数。 一个矩阵需要计 算（rows- 1 )次加窗互相关 （rows为矩阵行数）， 则由 （rows -1 )个 block实 现， 每个 block计算一个相邻两行的互相关系数。 因此内核函数线程数如下： 内核函数的线程数 ComputeXcross«< blocks, threads »>(), 其中， blocks = dim3(rows -1); threads = dim3(lagNum); lagNum是 lag (变化 间隔参数）总数。
假设一个窗口的互相关系数（有 lagNum个互相关系数）用 α」'表示， 通过 上述互相关计算后， 得到互相关系数矩阵包含（roWS-l)*WinNUm个 Cij , i = 0,1... rows -2 j = 0,1... WinNum- 1 步骤 206, 获取该互相关系数矩阵中最大互相关系数对应的变化间隔参 数值。
具体的， 采用内核函数 MatDelay对互相关系数矩阵求极值和位移。 求一 个 α」'内最大相关系数对应的 lag值 lag_
步骤 208, 根据该变化间隔参数值求取位移。
具体的， 由 lag_max求位移 delay delay的求取方法如公式（6):
delay = ( lag _ max + delta _ t ) * factor ( 6 ) 其中， delta_t采用公式（7)获取。 fact。r = c/(2*fe), c是声速 1540m/s, fs采样频率, *表示乘法。
当微小位移引起的超声回波信号时移 小于采样周期时， 可以通过互相 关系数求出， 公式如下：
St = 2K ， (7)
Zk n^ +\,s)-Zk n^ -\,s) 式 （7 ) 中， " 是最大互相关系数对应的坐标， Z是角度运算， 即 Z(a + ¾) = arctan(¾/a)? 得到了时移值 。 本实施例中， delta_t即为时移 可通 过式（7)计算得到。 通过上述计算得到位移矩阵。 输入数据矩阵大小为 rOWXCO/， 则得到的位 移矩阵为（row - 1) X WinNum , WinNum为窗口数。
此外 微小位移 d ,s、可通过超声波速 C求得： d(m, s) = St— ( 8 )
2
本 实 施 例 中 ， 内 核 函 数 线 程数 ： blocks=dim3(rows-l) , threads=dim3 (WinNum) , MatDelay «< blocks, threads »>(...), 其中 WinNum 窗个数， 通过这些线程一次计算后， 所有窗都计算了位移 delay。 互相关系数求位移可采用 CUDA的流方式求取。该 CUDA的流方式可以 隐藏数据传输时间， 进一步缩短程序运行时间。 例如在求位移的流处理示例 如下：
〃数据由内存输入显存， 这里是 4个矩阵， M=4, 这里代码并不完整。 cudaMemcpyAsync(d_datal ,h_datal , ... ,cudaMemcpyHostToDevice,stream[
0])；
cudaMemcpyAsync (d_data2 ,h_data2 , ...cudaMemcpyHostToDevice,strea m[l]);
cudaMemcpyAsync (d_data3 ,h_data3 , ...cudaMemcpyHostToDevice,strea m[2]);
cudaMemcpyAsync(d_data4,h_data4, ...cudaMemcpyHostToDevice,strea m[3]);
〃GPU计算互相关
ComputeXcross«<dBlock,dThread,0,stream[0]»>( datal , ...);
ComputeXcross«<dBlock,dThread,0,stream[l]»>( data2...);
ComputeXcross«<dBlock,dThread,0,stream[2]»>( data3...);
ComputeXcross«<dBlock,dThread,0,stream[3]»>( data4...);
//GPU计算移位
MatDelay«<blockForDelay,threadForDelay,0,stream[0]»>();
MatDelay«<blockForDelay,threadForDelay,0,stream[ 1 ]»>();
MatDelay«<blockForDelay,threadForDelay,0,stream[2]»>();
MatDelay«<blockForDelay,threadForDelay,0,stream[3]»>();
//结果从显存拷贝回内存 cudaMemcpyAsync (h data 1 ,
d datal ... ,cudaMemcpyDeviceToHost,stream[0]);
cudaMemcpyAsync (h_data2 ,
d_data2... ,cudaMemcpyDeviceToHost,stream[ 1 ]);
cudaMemcpyAsync (h_data3 ,
d_data3... ,cudaMemcpyDeviceToHost,stream[2]);
cudaMemcpyAsync(h_data4,
d_data4...,cudaMemcpyDeviceToHost,stream[3]);
II释放 stream
for (int i = 0 ; i < 4; ++i)
{
cudaStreamDestroy(stream[i]);
}
同样， 可采用该流的方式进行前向滤波处理、 逆向滤波处理和希尔伯特 变换。
此外， 还可对得到的预定数量位移矩阵进行平滑处理， 当信号噪音比较 大时， 可多次进行平滑处理。 步骤 108 , 位移合并步骤， 对每个位移矩形进行时间累加， 取累积和最 大的列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵。
具体的， 对得到的 m个位移矩阵分别进行处理： 进行时间 （列方向） 累 加， 取累积和最大的列位移， 再将 m个位移矩阵抽取的列位移合并， 生成新 的位移矩阵， 即 m X (row - 1)。
步骤 110, 插值步骤， 对得到的新的位移矩阵进行时间方向上的三次样 条插值。
该插值步骤包括： 采用中央处理器求取每个插值段之间的三次样条插值 函数的系数， 再采用图形处理器根据该三次样条插值函数并行计算每段的插 值倍数减 1个插值点的值， 其中， 图形处理器所调用的内核函数的 blocks数 量为待插值矩阵的行数， threads数量为待插值矩阵的列数。 三次样条插值倍 数可根据需要设定。 用三次样条插值方法在相邻两点间补充 p个点。 本实施 例中， 采用的插值倍数为 20 , 即相邻两点中间补充 19个点， p=19。
具体的， 插值段 [χ， χΜ ]的三次样条插值函数 Si (X)如下：
si (x) = f (x - xf )3 + (x _ x )2 + cf (x -xi ) r di ( 9 ) 式（9 ) 中， "'·、 b c 为插值函数系数， 可通过解线性方程获得， 因线性方程是不定解方程， 可以采用自然插值（Natural splines )确定方程解。
求插值系数的线性方程矩阵为三对角矩阵， 直接 LU分解求解。 由于 LU 直接分解是不断迭代过程， 采用 CPU计算出全部插值函数系数， 再用 GPU 并行计算出所有插值点的值。
设计内核函数 spline_kern计算一个 [χ^·"]插值段的所有值。 对输入数据 矩阵每行插值， 需要线程数如下： spline_kem«<blocks, threads»>(), 其中， blocks=dim3 (inputMat->rows)； threads=dim3 (inputMat->cols- 1 ) , inputMat->rows 是待插值矩阵行数， inputMat->cols是待插值矩阵的列数。
步骤 112, 拟合步骤， 对插值结果进行拉东变换， 得到拉东变换矩阵。 如图 3所示， 为拉动变换示意图。 图 3可知， 依次求一点到对端一点的累 加和， 考虑实际应用， 终点的序号不能小于起点的序号。 假设矩阵宽度为 Q, 需要计算 Q*(Q+l)/2次斜对角线的和。 为了存取数据方便， 通常 grid和 block的 维数直接用矩阵的长宽相关值， 这里经过插值后的输入矩阵， Q的值通常较 大比如 1960, 那么拉东变换矩阵大小为 1960*1960 , 而计算能力为 3.5的 GPU 支持 block内部最大线程数是 1024, 因此 block内部线程数不能用 Q定义。
定义 grid ( X*X ), block ( Y*Y )„ 其中， Υ取 16或者 32, Υ取满足 Χ*Υ>= Q最小值即可。 索引值按如下计算：
index—start = blockldx.x * blockDim.x + threadldx.x;
index—end = blockldx.y * blockDim.y + threadldx.y;
其中， *表示乘法。
即 grid的 x维和 block的 x维组合表示起始端点的序号 index_start ( 0到 1959 ), grid的 y维和 block的 y维组合表示终点端的序号 index_end。 内核函数实现一个斜对角求和（从初始点 index_start, 到终点 index_end 的对角连线上的数据求和）功能， 得到的结果 Sum ( index start, index end ) 保存到输出数据中。
步骤 114, 估算步骤， 获取该声辐射力脉冲成像的参数， 根据该参数及 拉东变换矩阵求取剪切波速率。
该估算步骤中该声辐射力脉冲成像的参数包括探头中两个阵元的间距、 探头个数、 输入数据矩阵的行数、 帧频和采样倍数。
该估算步骤中求取剪切波速率的公式为：
vs = pitch * pitchNum * (dataRow - \) I{L * fs * spline)
其中， 为剪切波速率， p fc/z是探头中两个阵元的间距， pitchNum 头个数， ^to tow是输入数据矩阵的行数， 是帧频， / ^是采样倍数， *表 示乘法， £ = (i?max).; - (i? max).x , 是位移， max是拉东变换矩阵最大值， L 是最大值对应的 j坐标减去最大值对应的 X坐标。
上述声辐射力脉冲成像估算方法， 通过对输入数据进行前向滤波处理和 逆向滤波处理， 消除了滤波造成的相对偏移， 并对滤波结果进行希尔伯特变 换， 减小了后续互相关计算位移时产生的误差， 使得计算结果更加准确， 且 对微小位移能精确估算， 提高了灵敏度。
另外， 估算过程中采用 CUDA并行计算实现， 提高了效率， 且在并行实 现上， 不增加重复计算量， 不使用线程同步和共享内存， 实现线程高度并行， 提高了处理效率。 如图 4所示， 为一个实施例中声辐射力脉冲成像估算系统的结构框图。 该声辐射力脉冲成像估算系统中未详尽描述的细节， 可参考上述对应的方法 实施例描述。 该声辐射力脉冲成像估算系统， 包括数据划分模块 410、 数据 滤波变换模块 420、 位移求取模块 430、 位移合并模块 440、 插值模块 450、 拟合模块 460和估算模块 470。 其中：
数据划分模块 410用于读取通过声辐射力脉冲成像所采集的输入数据， 并将该输入数据分成预定数量矩阵。 具体的， 将输入数据按侧向釆集位置分为 m个数据矩阵。 m为侧向釆集 位置的条数。每个数据矩阵大小为 rOW X CO/，表示包含 row帧数据（时间方向）， 每帧有 co/个元素 （深度方向）。
数据滤波变换模块 420用于分别对每个数据矩阵进行前向滤波处理和逆 向滤波处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据。
本实施例中， 1 )该数据滤波变换模块 420进行前向滤波处理的公式为： 当 i < iParalen时：
tOutPutik, ΐ) = ^ (param[j] * input k, i - j)) + ^ (param[j] * input(k,0))
=0 当 i > iParalen时：
tOutPutik, ΐ) = ^ (param[j] * input k, i - j)) 其中， param m 是滤波器长度， /¾rot为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， OWPW是滤波结果矩阵， *表示乘法。
该前向滤波处理釆用 GPU 并行处理。 该前向滤波处理的内核函数 Bandpass front (即图形处理器调用的内核函数） 调用的线程数设计： Bandpass_front«< blocks, threads»>( input, aram, iParaLen, tOutPut) , 此 处 blocks = dim3(input->rows); threads = dim3(input->cols); 其中, input ->rows 是输入数据矩阵的行数， input->cols 是输入数据矩阵的列数。 该输入数据矩 阵即为待滤波矩阵。
2 )该数据滤波变换模块 420进行逆向滤波处理的公式为：
当 i > iParalen时：
output (Ji, Width - l - i) = ^ (param[j] * input k, Width _ 1 + j - i))
=o 其中， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， 。i ^^为逆向滤波的滤波结果矩阵， 为滤波器宽度， *表示乘法。
该逆向滤波处理釆用图形处理器处理。 该逆向滤波处理的内核函数
Bandpass back (即图形处理器调用的内核函数） 调用的线程数设计： Bandpass_back«< blocks , threads»>() , 其中, blocks = dim3 (input ->rows); threads = dim3 (input->cols) , 滤波器和前向滤波器相同， input->rows是输入数 据矩阵的行数， input->cols是输入数据矩阵的列数。 该输入数据矩阵即为待滤 波矩阵。 内核函数 Bandpass_back将一个元素按照公式 （3 ) 处理， 将获取的 结果存放在 outPut对应的位置。
3 )该数据滤波变换模块 420进行希尔伯特变换公式为：
rA (m, t) = r(m, t) * (l /(7rt)) ( 4 ) 式（4 ) 中， Λ("^)为信号 r(w 的希尔伯特变换， w为采样点的坐标， 是时间， *是卷积运算。
在算法上采用傅里叶变换和逆变换来实现希尔伯特变换， 具体为： 对初 始的矩阵进行傅里叶变换， 将变换结果的实部乘 -1后和虚部交换， 再进行傅 里叶归一化逆变换， 再交换实虚部， 再和初始的矩阵相加得到解析数据。
希尔伯特变换采用 GPU优化设计， 其中， 傅里叶变换和逆变换均使用 CUDA (即 CPU+GPU编程模型 ) 的 CUFFT库函数 cufftExecC2C。
其中， 两次实虚部交换以及矩阵相加用 2个内核函数来实现， 第一个内核 函数 kerne 1_ 1实现第一次的傅里叶变换前的实部虚部交换， 第二个内核函数 kernel_2实现第二次实部虚部交换， 并且交换后相加， 因 cufftExecC2C结果没 有归一化处理， 在第二次实虚部交换时候处理。
内核函数调用顺序： 1 ) cufftExecC2C (正变换）； 2 ) kemel l； 3 ) cufftExecC2C (逆变换）； 4 ) kernel_2。
两次内核函数的 blocks数量和 threads数量相同，均为 blocks = dim3 (input -〉 rows); threads = dim3 (input -〉 cols), input -〉 rows是输入数据矩阵的行数, input ->cols是输入数据矩阵的列数。
在一个实施例中， 该数据滤波变换模块 420还采用多个流按顺序依次对 预定数量数据矩阵进行前向滤波处理、 逆向滤波处理和希尔伯特变换处理。 具体的， 因 GPU计算的一个特点是待处理数据需要从内存传送到显存， 结果 数据需要从显存传回内存这两个数据传输过程。 由 CPU调用 GPU函数接口 把待处理数据传输到显存， GPU通过显存读取数据， 计算数据， 把结果存放 显存， CPU再调用数据传输函数接口将显存数据传回内存， 进行下一步逻辑 处理。 数据在 CPU和 GPU之间的传输过程可与 GPU计算过程并行。 可使用 多个流（stream ) 实现数据传输与 GPU计算过程并行处理。 流与流之间是并 行运行的， 因此虽然 GPU—次只能处理一个内核函数， 而当多个流存在时， GPU在处理内核函数的同时， 显存可以传输数据， 从而提高设备的利用率。
位移求取模块 430用于对该解析数据进行互相关运算， 求取得到预定数 量位移矩阵。
如图 5所示， 该位移求取模块 430包括选取子模块 432、 互相关系数计 算子模块 434、 比较子模块 436和位移计算子模块 438。 其中：
选取子模块 432用于选取窗口， 该窗口长度为预定值。 本实施例中， 窗 口长度为 M , M为正数。
互相关系数计算子模块 434用于根据该窗口长度按步长滑动的方式对该 解析数据计算每个窗口的互相关系数， 并形成互相关系数矩阵。
本实施例中， 采用加窗互相关求取位移， 加窗互相关是指一次只针对长 度为窗口长度的一段求互相关， 窗口位置随步长滑动， 直到数据末位为止。
每个窗口的互相关系数计算公式为：
k(q, s) = ^ rA (m + u, s)rA (m + u + q, s + \) ( 5 ) 式（5 )中， 是互相关系数， + M， 是解析数据， r m + u + q, s + \) 是解析数据 + Μ + ^ + 求复共轭， s是采样时间， g为 [-Μ,Μ]。
g每一次变化称一个 lag (即变化间隔参数）， lag的范围 [-Μ,Μ]。 当估算 位移比较小时， lag取一个范围 [-N, N] , N < M。
采用内核函数 ComputeXcross计算所有窗（ WinNum个 )中一次 lag的互 相关系数。 WinNum=(InputWidth - WindowHW)/Step+l , 其中， WindowHW 是窗口宽度， 步长为 Step, InputWidth为输入数据矩阵的列数， 为了减少重 复计算提高时间效率， WinNum个互相关系数计算方法如下： 计算数据矩阵 相邻两行 Row(i), Row(i+l)互相关系数， 则先计算第 0个窗口的互相关系数。 第 n个窗口的互相关系数计算为：第 n-1个窗口的互相关系数减去第 n-1个窗 口前 Step个元素的互相关系数，加上第 n个窗口最后 Step个元素的互相关系 数。 例如， 第一个窗口的互相关系数为： 第 0个窗口的互相关系数减去第 0 个窗口前 Step个元素的互相关系数，加上第一个窗口最后 Step个元素的互相 关系数； 第二个窗口的互相关系数： 第一个窗口的互相关系数减去第一个窗 前 Step个元素的互相关系数，加上第二个窗口最后 Step个元素的互相关系数， 依次计算到最后一个窗口为止。
内核函数 ComputeXcross计算 1个 lag互相关系数， 设计 1个 block的线程数 等于 lag总数， 可以计算出相邻两行所有窗口的互相关系数。 个矩阵需要计 算（rows- 1 )次加窗互相关 （rows为矩阵行数）， 则由 （rows -1 )个 block实 现， 每个 block计算一个相邻两行的互相关系数。 因此内核函数线程数如下： 内核函数的线程数 ComputeXcross«< blocks, threads »>(), 其中， blocks = dim3(rows -1); threads = dim3(lagNum); lagNum是 lag (变化 间隔参数）总数。
假设一个窗口的互相关系数（有 lagNum个互相关系数）用 α」'表示， 通过 上述互相关计算后， 得到互相关系数矩阵包含（r0WS-l)*WinNUm个 Cij , i = 0,1... rows -2 j = 0,1... WinNum- 1 比较子模块 436用于获取该互相关系数矩阵中最大互相关系数对应的变 化间隔参数值。 具体的， 采用内核函数 MatDelay对互相关系数矩阵求极值和 位移。 求一个 α」'内最大相关系数对应的 lag值 lag_max
位移计算子模块 438用于根据该变化间隔参数值求取位移。
具体的， 由 lag_max求位移 delay delay的求取方法如公式（6):
delay = ( lag _ max + delta _ t ) * factor ( 6 ) 其中， delta_t采用公式（7)获取。 &ct。r = c/(2*fe), c是声速 1540m/s, fs采样频率, *表示乘法。
当微小位移引起的超声回波信号时移 小于采样周期时， 可以通过互相 关系数求出， 公式如下：
8t = 2K ， (7)
Zk nm^ +\,s)-Zk nm^ -\,s) 式 （7 ) 中， " 是最大互相关系数对应的坐标， Z是角度运算， 即 Z(a + ¾) = arctan(¾/a)? 得到了时移值 。 本实施例中， delta t即为时移 可通 过式（7 )计算得到。 通过上述计算得到位移矩阵。 输入数据矩阵大小为 rOWX CO/， 则得到的位 移矩阵为（row - 1) X WinNum , WinNum为窗口数。
此外 微小位移 d(m,S)可通过超声波速 C求得：
d(m,s) = St— ( 8 )
2
本 实 施 例 中 ， 内 核 函 数 线 程数 ： blocks=dim3(rows-l) , threads=dim3 (WinNum) , MatDelay «< blocks, threads >>>(), 其中 WinNum 窗个数， 通过这些线程一次计算后， 所有窗都计算了位移 delay 互相关系数求位移可采用 CUDA的流方式求取。
位移合并模块 440用于对每个位移矩形进行时间累加， 取累积和最大的 列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵。 具 体的， 对得到的 m个位移矩阵分别进行处理： 进行时间 （列方向） 累加， 取 累积和最大的列位移， 再将 m个位移矩阵抽取的列位移合并， 生成新的位移 矩阵， 即 mx 1)。
插值模块 450用于对得到的新的位移矩阵进行时间方向上的三次样条插 值。
本实施例中， 该插值模块 450还用于采用中央处理器求取每个插值段之 间的三次样条插值函数的系数， 再采用图形处理器根据该三次样条插值函数 并行计算每段的插值倍数减 1个插值点的值， 其中， 图形处理器所调用的内 核函数的 blocks数量为待插值矩阵的行数, threads数量为待插值矩阵的列数。 三次样条插值倍数可根据需要设定。 用三次样条插值方法在相邻两点间补充 p个点。 本实施例中， 采用的插值倍数为 20, 即相邻两点中间补充 19个点， p=19
具体的 插值段 [χ χΜ ]的三次样条插值函数 Si (X)如下：
si (x) = f(x-xf)3 + (x_x )2 +cf(x -xi) rdi ( 9 ) 式（9) 中， a'、 b 为插值函数系数， 可通过解线性方程获得， 因线性方程是不定解方程， 可以釆用自然插值（Natural splines )确定方程解。 求插值系数的线性方程矩阵为三对角矩阵， 直接 LU分解求解。 由于 LU 直接分解是不断迭代过程， 釆用 CPU计算出全部插值函数系数， 再用 GPU 并行计算出所有插值点的值。
设计内核函数 spline_kern计算一个 [χ^·"]插值段的所有值。 对输入数据 矩阵每行插值， 需要线程数如下： spline_kem«<blocks, threads»>(), 其中， blocks=dim3 (inputMat->rows) ； threads=dim3 (inputMat->cols -1) , inputMat->rows是待插值矩阵行数， inputMat->cols是待插值矩阵的列数。
拟合模块 460用于对插值结果进行拉东变换， 得到拉东变换矩阵。
如图 3所示， 为拉动变换示意图。 图 3可知， 依次求一点到对端一点的累 加和， 考虑实际应用， 终点的序号不能小于起点的序号。 假设矩阵宽度为 Q, 需要计算 Q*(Q+l)/2次斜对角线的和。 为了存取数据方便， 通常 grid和 block的 维数直接用矩阵的长宽相关值， 这里经过插值后的输入矩阵， Q的值通常较 大比如 1960, 那么拉东变换矩阵大小为 1960*1960 , 而计算能力为 3.5的 GPU 支持 block内部最大线程数是 1024, 因此 block内部线程数不能用 Q定义。
定义 grid ( X*X ), block ( Y*Y )„ 其中， Υ取 16或者 32, Υ取满足 Χ*Υ>= Q最小值即可。 索引值按如下计算：
index—start = blockldx.x * blockDim.x + threadldx.x;
index—end = blockldx.y * blockDim.y + threadldx.y;
其中， *表示乘法。 即 grid的 x维和 block的 x维组合表示起始端点的序号 index start ( 0到 1959 ) , grid的 y维和 block的 y维组合表示终点端的序号 index—end。
内核函数实现一个斜对角求和（从初始点 index_start, 到终点 index_end 的对角连线上的数据求和）功能， 得到的结果 Sum ( index start, index end ) 保存到输出数据中。
估算模块 470用于获取该声辐射力脉冲成像的参数， 根据该参数及拉东 变换矩阵求取剪切波速率。
本实施例中， 该声辐射力脉冲成像的参数包括探头中两个阵元的间距、 探头个数、 输入数据矩阵的行数、 帧频和釆样倍数；
该估算模块用于求取剪切波速率的公式为：
vs = pitch * pitchNum * (dataRow - \) I{L * fs * spline)
其中， 为剪切波速率， p fc/z是探头中两个阵元的间距， pitchNum 头个数， ^to tow是输入数据矩阵的行数， 是帧频， / ^是釆样倍数， *表 示乘法， £ = (i?max).; - (i? max).x , 是位移， max是拉东变换矩阵最大值， L 是最大值对应的 j坐标减去最大值对应的 X坐标。
上述声辐射力脉冲成像估算系统， 通过对输入数据进行前向滤波处理和 逆向滤波处理， 消除了滤波造成的相对偏移， 并对滤波结果进行希尔伯特变 换， 减小了后续互相关计算位移时产生的误差， 使得计算结果更加准确， 且 对微小位移能精确估算， 提高了灵敏度。
另外， 估算过程中釆用 CUDA并行计算实现， 提高了效率， 且在并行实 现上， 不增加重复计算量， 不使用线程同步和共享内存， 实现线程高度并行， 提高了处理效率。
上述声辐射力脉冲成像估算方法和系统， 模拟结果： 实验中给定仿体的 剪切波速度为 3m/s (米 /秒）， 计算结果为 3.002m/ s, 误差仅为 0.002m/s。
GPU优化效果： 本估算方法以 4组 100*512的釆集数据作为一组输入数 据， 设置窗口长度 100, lag总数为 81 , —次计算耗时为 80ms。 该耗时时间 包括了 GPU与 CPU之间数据传输时间和 GPU计算时间， 不包括 GPU申请 显存时间。 由于 GPU申请显存比较慢（申请足够的内存用时大约 500ms ), 考虑到申请的显存可以重复使用， 可在系统启动时一次性申请足够显存， 在 系统关闭时释放； 申请内存时间不作为运行时间计算是合理的。对比 GPU优 化前， 利用配置为双核主频 3.3G、 3.29G内存 2G的个人计算机运行算法， 耗 时为 2000ms左右。
请参考图 6, 为本发明的声辐射力脉冲成像估算方法和系统所运行的终 端的结构示意图， 该终端可以用于实施上述实施例中提供的声辐射力脉冲成 像估算方法。 具体来讲：
终端 1200 可以包括有一个或一个以上计算机可读存储介质的存储器 110、 输入单元 120、 显示单元 130、 音频电路 140、 包括有一个或者一个以 上处理核心的处理器 150 (处理器 140包括中央处理器（CPU )和图形处理器 ( GPU ) )、 以及电源 160等部件。 本领域技术人员可以理解， 图 6中示出的 终端结构并不构成对终端的限定， 可以包括比图示更多或更少的部件， 或者 组合某些部件， 或者不同的部件布置。 其中：
存储器 110可用于存储软件程序以及模块， 处理器 150通过运行存储在 存储器 110的软件程序以及模块， 从而执行各种功能应用以及数据处理。 存 储器 110可主要包括存储程序区和存储数据区， 其中， 存储程序区可存储操 作系统、 至少一个功能所需的应用程序 （比如声音播放功能、 图像播放功能 等）等； 存储数据区可存储根据终端 1200的使用所创建的数据（比如声辐射 力脉冲成像数据）等。 此外， 存储器 110可以包括高速随机存取存储器， 还 可以包括非易失性存储器， 例如至少一个磁盘存储器件、 闪存器件、 或其他 易失性固态存储器件。 相应地， 存储器 110还可以包括存储器控制器， 以提 供处理器 150和输入单元 120对存储器 110的访问。
输入单元 120可用于接收输入的数字或字符信息， 以及产生与用户设置 以及功能控制有关的键盘、 鼠标、 操作杆、 光学或者轨迹球信号输入。 具体 地，输入单元 120可包括触敏表面 121以及其他输入设备 122。触敏表面 121 , 也称为触摸显示屏或者触控板， 可收集用户在其上或附近的触摸操作 （比如 用户使用手指、 触笔等任何适合的物体或附件在触敏表面 121上或在触敏表 面 121附近的操作）， 并根据预先设定的程式驱动相应的连接装置。 可选的， 触敏表面 121可包括触摸检测装置和触摸控制器两个部分。 其中， 触摸检测 装置检测用户的触摸方位， 并检测触摸操作带来的信号， 将信号传送给触摸 控制器； 触摸控制器从触摸检测装置上接收触摸信息， 并将它转换成触点坐 标， 再送给处理器 150, 并能接收处理器 150发来的命令并加以执行。 此外， 可以采用电阻式、电容式、红外线以及表面声波等多种类型实现触敏表面 121。 除了触敏表面 121 , 输入单元 120还可以包括其他输入设备 122。 具体地， 其 他输入设备 122可以包括但不限于物理键盘、 功能键（比如音量控制按键、 开关按键等）、 轨迹球、 鼠标、 操作杆等中的一种或多种。
显示单元 130可用于显示由用户输入的信息或提供给用户的信息以及终 端 1200的各种图形用户接口， 这些图形用户接口可以由图形、 文本、 图标、 视频和其任意组合来构成。 显示单元 130可包括显示面板 131 , 可选的， 可 以采用 LCD(Liquid Crystal Display , 液晶显示器）、 OLED(Organic Light-Emitting Diode,有机发光二极管)等形式来配置显示面板 131。进一步的， 触敏表面 121可覆盖显示面板 131 , 当触敏表面 121检测到在其上或附近的 触摸操作后， 传送给处理器 150以确定触摸事件的类型， 随后处理器 150根 据触摸事件的类型在显示面板 131上提供相应的视觉输出。 虽然在图 6中， 触敏表面 121与显示面板 131是作为两个独立的部件来实现输入和输入功能， 但是在某些实施例中， 可以将触敏表面 121与显示面板 131集成而实现输入 和输出功能。
音频电路 140、 扬声器 141 , 传声器 142可提供用户与终端 1200之间的 音频接口。 音频电路 140可将接收到的音频数据转换后的电信号， 传输到扬 声器 141 , 由扬声器 141转换为声音信号输出； 另一方面， 传声器 142将收 集的声音信号转换为电信号， 由音频电路 140接收后转换为音频数据， 再将 音频数据输出处理器 150处理。 音频电路 140还可能包括耳塞插孔， 以提供 外设耳机与终端 1200的通信。
处理器 150是终端 1200的控制中心，利用各种接口和线路连接整个计算 机的各个部分， 通过运行或执行存储在存储器 110内的软件程序和 /或模块， 以及调用存储在存储器 110内的数据，执行终端 1200的各种功能和处理数据。 可选的， 处理器 150可包括一个或多个处理核心； 优选的， 处理器 150可集 成应用处理器和调制解调处理器， 其中， 应用处理器主要处理操作系统、 用 户界面和应用程序等， 调制解调处理器主要处理无线通信。 可以理解的是， 上述调制解调处理器也可以不集成到处理器 150中。
终端 1200还包括给各个部件供电的电源 160 (比如电池）， 优选的， 电 源可以通过电源管理系统与处理器 150逻辑相连， 从而通过电源管理系统实 现管理充电、 放电、 以及功耗管理等功能。 电源 160还可以包括一个或一个 以上的直流或交流电源、 再充电系统、 电源故障检测电路、 电源转换器或者 逆变器、 电源状态指示器等任意组件。
尽管未示出，终端 1200还可以包括摄像头、蓝牙模块等，在此不再赘述。 该终端 1200为计算机。
本发明还提供了一个或多个包含计算机可执行指令的计算机存储介质， 该计算机可执行指令包含一个或一个以上程序， 所述一个或一个以上程序被 一个或一个以上的终端的处理器用来执行一种声辐射力脉冲成像估算方法， 所述方法包括以下步骤：
数据划分步骤， 读取通过声辐射力脉冲成像所釆集的输入数据， 并将所 述输入数据分成预定数量数据矩阵；
数据滤波变换步骤， 分别对每个数据矩阵进行前向滤波处理和逆向滤波 处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；
位移求取步骤， 对所述解析数据进行互相关运算， 求取得到预定数量位 移矩阵；
位移合并步骤， 对每个位移矩形进行时间累加， 取累积和最大的列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值步骤， 对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合步骤， 对插值结果进行拉东变换， 得到拉东变换矩阵；
估算步骤， 获取所述声辐射力脉冲成像的参数， 根据所述参数及拉东变 换矩阵求取剪切波速率。
在其中一个实施例中， 所述数据滤波变换步骤中， 釆用多个流按顺序依 次对预定数量数据矩阵进行前向滤波处理、 逆向滤波处理和希尔伯特变换处 理。
在其中一个实施例中，所述数据滤波变换步骤中前向滤波处理的公式为： 当 i < iParalen时：
tOutPut{k, ί) - ^ (param[j] * input{k, i j)) + ^ (param[j] * input kfi)) 当 〉 iParalen时：
tOutPut{k, i) - ^ (param[j] * input{k, i j)) 其中， ra 是滤波器， 是滤波器长度， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， C ^是滤波结果矩阵， *表示乘法； 所述数据滤波变换步骤中逆向滤波处理的公式为：
当 i > iParalen时：
output ( i, Width - 1 - /) = ^ (param[j] * input{k, Width _ 1 + j— /)) 其中， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， 。i ^^为逆向滤波的滤波结果矩阵， 为滤波器宽度， *表示乘法；
所述前向滤波处理的内核函数和逆向滤波处理均采用图形处理器并行处 理， 所述图形处理器调用的内核函数的 blocks数量和 threads数量相同， 所述 blocks数量为待滤波矩阵的行数， 所述 threads数量为待滤波矩阵的列数。
在其中一个实施例中， 所述数据滤波变换步骤中希尔伯特变换包括： 对初始的数据矩阵进行傅里叶变换， 将变换结果的实部乘 -1 后和虚部交 换， 再进行傅里叶归一化逆变换， 再交换实虚部， 再和初始的数据矩阵相加 得到解析数据。
在其中一个实施例中， 所述位移求取步骤包括：
选取窗口， 所述窗口长度为预定值；
根据所述窗口长度按步长滑动的方式对所述解析数据计算每个窗口的互 相关系数， 并形成互相关系数矩阵；
获取所述互相关系数矩阵中最大互相关系数对应的变化间隔参数值； 根据所述变化间隔参数值求取位移。
在其中一个实施例中， 所述插值步骤包括：
采用中央处理器求取每个插值段之间的三次样条插值函数的系数， 再采 用图形处理器根据所述三次样条插值函数并行计算每段的插值倍数减 1个插 值点的值， 其中， 图形处理器所调用的内核函数的 blocks数量为待插值矩阵 的行数， threads数量为待插值矩阵的列数。
在其中一个实施例中， 所述估算步骤中所述声辐射力脉冲成像的参数包 括探头中两个阵元的间距、 探头个数、 输入数据矩阵的行数、 帧频和采样倍 数；
所述估算步骤中求取剪切波速率的公式为：
vs = pitch * pitchNum * (dataRow - \) I{L * fs * spline)
其中， 为剪切波速率， p fc/z是探头中两个阵元的间距， pitchNum 头个数， ^to tow是输入数据矩阵的行数， 是帧频， / ^是采样倍数， *表 示乘法， £ = (i?max).; - (i? max).x , 是位移， max是拉东变换矩阵最大值， L 是最大值对应的 y坐标减去最大值对应的 X坐标。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流 程， 是可以通过计算机程序来指令相关的硬件来完成， 所述的程序可存储于 一计算机可读取存储介质中， 该程序在执行时， 可包括如上述各方法的实施 例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（ Read-Only Memory, ROM )或随机存储记忆体 ( Random Access Memory, RAM )等。 细， 但并不能因此而理解为对本发明专利范围的限制。 应当指出的是， 对于 本领域的普通技术人员来说， 在不脱离本发明构思的前提下， 还可以做出若 干变形和改进， 这些都属于本发明的保护范围。 因此， 本发明专利的保护范 围应以所附权利要求为准。

Claims

权利要求书
1、 一种声辐射力脉冲成像估算方法， 运行于计算机上， 其特征在于， 包 括：
数据划分步骤， 读取通过声辐射力脉冲成像所采集的输入数据， 并将所 述输入数据分成预定数量数据矩阵；
数据滤波变换步骤， 分别对每个数据矩阵进行前向滤波处理和逆向滤波 处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；
位移求取步骤， 对所述解析数据进行互相关运算， 求取得到预定数量位 移矩阵；
位移合并步骤， 对每个位移矩形进行时间累加， 取累积和最大的列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值步骤， 对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合步骤， 对插值结果进行拉东变换， 得到拉东变换矩阵；
估算步骤， 获取所述声辐射力脉冲成像的参数， 根据所述参数及拉东变 换矩阵求取剪切波速率。
2、 根据权利要求 1所述的声辐射力脉冲成像估算方法， 其特征在于， 所 述数据滤波变换步骤中， 采用多个流按顺序依次对预定数量数据矩阵进行前 向滤波处理、 逆向滤波处理和希尔伯特变换处理。
3、 根据权利要求 1所述的声辐射力脉冲成像估算方法， 其特征在于， 所 述数据滤波变换步骤中前向滤波处理的公式为：
当 i < iParalen时：
tOutPut{k, ί) - ^ (param[j] * input{k, i j)) + ^ (param[j] * input kfi)) 当 〉 iParalen时：
tOutPut{k, i) - ^ (param[j] * input{k, i j)) 其中， ra 是滤波器， 是滤波器长度， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， C ^是滤波结果矩阵， *表示乘法； 所述数据滤波变换步骤中逆向滤波处理的公式为： 当 i > iParalen时： output (Ji, Width - l - i) = ^ (param[j] * input k, Width _ 1 + j - i))
=o 其中， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， 。i ^^为逆向滤波的滤波结果矩阵， 为滤波器宽度， *表示乘法；
所述前向滤波处理的内核函数和逆向滤波处理均采用图形处理器并行处 理， 所述图形处理器调用的内核函数的 blocks数量和 threads数量相同， 所述 blocks数量为待滤波矩阵的行数， 所述 threads数量为待滤波矩阵的列数。
4、 根据权利要求 1所述的声辐射力脉冲成像估算方法， 其特征在于， 所 述数据滤波变换步骤中希尔伯特变换包括：
对初始的数据矩阵进行傅里叶变换， 将变换结果的实部乘 -1 后和虚部交 换， 再进行傅里叶归一化逆变换， 再交换实虚部， 再和初始的数据矩阵相加 得到解析数据。
5、 根据权利要求 1所述的声辐射力脉冲成像估算方法， 其特征在于， 所 述位移求取步骤包括：
选取窗口， 所述窗口长度为预定值；
根据所述窗口长度按步长滑动的方式对所述解析数据计算每个窗口的互 相关系数， 并形成互相关系数矩阵；
获取所述互相关系数矩阵中最大互相关系数对应的变化间隔参数值； 根据所述变化间隔参数值求取位移。
6、 根据权利要求 1所述的声辐射力脉冲成像估算方法， 其特征在于， 所 述插值步骤包括：
采用中央处理器求取每个插值段之间的三次样条插值函数的系数， 再采 用图形处理器根据所述三次样条插值函数并行计算每段的插值倍数减 1个插 值点的值， 其中， 图形处理器所调用的内核函数的 blocks数量为待插值矩阵 的行数， threads数量为待插值矩阵的列数。
7、 根据权利要求 1所述的声辐射力脉冲成像估算方法， 其特征在于， 所 述估算步骤中所述声辐射力脉冲成像的参数包括探头中两个阵元的间距、 探 头个数、 输入数据矩阵的行数、 帧频和采样倍数； 所述估算步骤中求取剪切波速率的公式为：
vs = pitch * pitchNum * (dataRow - \) I{L * fs * spline)
其中， 为剪切波速率， p fc/z是探头中两个阵元的间距， pitchNum 个数， ^to tow是输入数据矩阵的行数， 是帧频， / ^是采样倍数， *表示 乘法， £ = (i?max).;— (i? max).x , 是位移， ? max是拉东变换矩阵最大值， 是 最大值对应的 j坐标减去最大值对应的 X坐标。
8、 一种声辐射力脉冲成像估算系统， 运行于计算机上， 其特征在于， 包 括：
数据划分模块， 用于读取通过声辐射力脉冲成像所采集的输入数据， 并 将所述输入数据分成预定数量矩阵；
数据滤波变换模块， 用于分别对每个数据矩阵进行前向滤波处理和逆向 滤波处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据； 位移求取模块， 用于对所述解析数据进行互相关运算， 求取得到预定数 量位移矩阵；
位移合并模块， 用于对每个位移矩形进行时间累加， 取累积和最大的列 位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值模块，用于对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合模块， 用于对插值结果进行拉东变换， 得到拉东变换矩阵； 估算模块， 用于获取所述声辐射力脉冲成像的参数， 根据所述参数及拉 东变换矩阵求取剪切波速率。
9、 根据权利要求 8所述的声辐射力脉冲成像估算系统， 其特征在于， 所 述数据滤波变换模块还用于采用多个流按顺序依次对预定数量数据矩阵进行 前向滤波处理、 逆向滤波处理和希尔伯特变换处理。
10、 根据权利要求 8所述的声辐射力脉冲成像估算系统， 其特征在于， 所述数据滤波变换模块进行前向滤波处理的公式为：
当 i < iParalen时：
2
tOutPut{k, ί) - ^ (param[j] * input{k, i j)) + ^ (param[j] * input kfi)) 当 〉 iParalen时： tOutPutik, ΐ) = ^ (param[j] * input k, i - j)) 其中， param m 是滤波器长度， /¾rot为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， OWPW是滤波结果矩阵， *表示乘法； 所述数据滤波变换模块进行逆向滤波处理的公式为：
当 i > iParalen时：
output (Ji, Width - l - i) = ^ (param[j] * input k, Width _ 1 + j - i))
=o 其中， 为待滤波矩阵， k、 分别表示待滤波矩阵的第 行和第 列， 。i ^^为逆向滤波的滤波结果矩阵， 为滤波器宽度， *表示乘法；
所述前向滤波处理的内核函数和逆向滤波处理均采用图形处理器并行处 理， 所述图形处理器调用的内核函数的 blocks数量和 threads数量相同， 所述 blocks数量为待滤波矩阵的行数， 所述 threads数量为待滤波矩阵的列数。
11、 根据权利要求 8所述的声辐射力脉冲成像估算系统， 其特征在于， 所述数据滤波变换模块进行希尔伯特变换包括：
对初始的数据矩阵进行傅里叶变换， 将变换结果的实部乘 -1 后和虚部交 换， 再进行傅里叶归一化逆变换， 再交换实虚部， 再和初始的数据矩阵相加 得到解析数据。
12、 根据权利要求 8所述的声辐射力脉冲成像估算系统， 其特征在于， 所述位移求取模块包括：
选取子模块， 用于选取窗口， 所述窗口长度为预定值；
互相关系数计算子模块， 用于根据所述窗口长度按步长滑动的方式对所 述解析数据计算每个窗口的互相关系数， 并形成互相关系数矩阵；
比较子模块， 用于获取所述互相关系数矩阵中最大互相关系数对应的变 化间隔参数值；
位移计算子模块， 用于根据所述变化间隔参数值求取位移。
13、 根据权利要求 8所述的声辐射力脉冲成像估算系统， 其特征在于， 所述插值模块还用于采用中央处理器求取每个插值段之间的三次样条插值函 数的系数， 再采用图形处理器根据所述三次样条插值函数并行计算每段的插 值倍数减 1个插值点的值， 其中， 图形处理器所调用的内核函数的 block数 为待插值矩阵的行数， thread数为待插值矩阵的列数。
14、 根据权利要求 8所述的声辐射力脉冲成像估算系统， 其特征在于， 所述声辐射力脉冲成像的参数包括探头中两个阵元的间距、 探头个数、 输入 数据矩阵的行数、 帧频和采样倍数；
所述估算模块用于求取剪切波速率的公式为：
vs = pitch * pitchNum * (dataRow - \) I{L * fs * spline)
其中， 为剪切波速率， p fc/z是探头中两个阵元的间距， pitchNum 头个数， ^to tow是输入数据矩阵的行数， 是帧频， / ^是采样倍数， *表 示乘法， £ = (i?max).; - (i? max).x , 是位移， max是拉东变换矩阵最大值， L 是最大值对应的 j坐标减去最大值对应的 X坐标。
15、 一个或多个包含计算机可执行指令的计算机存储介质， 所述计算机 可执行指令用于执行一种声辐射力脉冲成像估算方法， 其特征在于， 所述方 法包括：
数据划分步骤， 读取通过声辐射力脉冲成像所采集的输入数据， 并将所 述输入数据分成预定数量数据矩阵；
数据滤波变换步骤， 分别对每个数据矩阵进行前向滤波处理和逆向滤波 处理， 并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；
位移求取步骤， 对所述解析数据进行互相关运算， 求取得到预定数量位 移矩阵；
位移合并步骤， 对每个位移矩形进行时间累加， 取累积和最大的列位移， 再将预定数量位移矩阵抽取的列位移合并， 生成新的位移矩阵；
插值步骤， 对得到的新的位移矩阵进行时间方向上的三次样条插值； 拟合步骤， 对插值结果进行拉东变换， 得到拉东变换矩阵；
估算步骤， 获取所述声辐射力脉冲成像的参数， 根据所述参数及拉东变 换矩阵求取剪切波速率。
PCT/CN2013/083049 2013-09-06 2013-09-06 声辐射力脉冲成像估算方法和系统、计算机存储介质 Ceased WO2015032064A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/CN2013/083049 WO2015032064A1 (zh) 2013-09-06 2013-09-06 声辐射力脉冲成像估算方法和系统、计算机存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2013/083049 WO2015032064A1 (zh) 2013-09-06 2013-09-06 声辐射力脉冲成像估算方法和系统、计算机存储介质

Publications (1)

Publication Number Publication Date
WO2015032064A1 true WO2015032064A1 (zh) 2015-03-12

Family

ID=52627706

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2013/083049 Ceased WO2015032064A1 (zh) 2013-09-06 2013-09-06 声辐射力脉冲成像估算方法和系统、计算机存储介质

Country Status (1)

Country Link
WO (1) WO2015032064A1 (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103431874A (zh) * 2013-09-06 2013-12-11 中国科学院深圳先进技术研究院 声辐射力脉冲成像估算方法和系统

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103431874A (zh) * 2013-09-06 2013-12-11 中国科学院深圳先进技术研究院 声辐射力脉冲成像估算方法和系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
QIU, GUANG: "Ultrasonic-based Tissue Elasticity Imaging Algorithm Study", CHINA MASTER'S THESES FULL-TEXT DATABASE, NO. 1, MEDICINE & PUBLIC HEALTH, no. 1, 15 January 2013 (2013-01-15), pages E060 - 18 *
ROSENZWEIG, S . ET AL.: "GPU-based Real-time Small Displacement Estimation with Ultrasound", IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, vol. 58, no. 2, 22 February 2011 (2011-02-22), pages 399 - 405 *
XING, YUANDING ET AL.: "Application of Zero-phase Filter in Dynamic Displacement Signal Processing Based on Laser Doppler Effect", ELECTRONIC MEASUREMENT TECHNOLOGY, vol. 32, no. 6, 30 June 2009 (2009-06-30), pages 48 - 50 *

Similar Documents

Publication Publication Date Title
CN103431874B (zh) 声辐射力脉冲成像估算方法和系统
CN102460356A (zh) 用于确定包括两个或更多个换能器的物体的表面上的一个或多个触击或触摸的位置的方法
CN109979469B (zh) 信号处理方法、设备及存储介质
EP3563909B1 (en) Step counting method and device for treadmill
Cuomo et al. A novel O (n) numerical scheme for ECG signal denoising
CN111383653A (zh) 语音处理方法及装置、存储介质、机器人
CN115795282A (zh) 激波管动态压力重构方法、装置、电子设备及存储介质
TW201216082A (en) Apparatus for adaptive time-frequency signal anaysis and method therefor
CN109186752B (zh) 基于图形处理器的水下声学信号采集、传输和检测系统
CN110755055A (zh) 一种用于确定脉搏波形的波形评价信息的方法与设备
CN114757229A (zh) 信号处理方法、装置、电子设备及介质
CN114611349B (zh) 一种基于隐式离散模型的频率域声波仿真方法及系统
CN113406610A (zh) 一种目标检测方法、装置、设备和存储介质
WO2015032064A1 (zh) 声辐射力脉冲成像估算方法和系统、计算机存储介质
Kajihara et al. Hardware implementation of autoregressive model estimation using Burg’s method for low-energy spectral analysis
WO2023216058A1 (zh) 一种信号起始点检测方法、装置、存储介质及电子设备
Ahmed et al. Towards smart health monitoring system for elderly people
CN110852307B (zh) 基于脑电信号的脑机接口的检测方法、系统、介质及设备
CN108742581B (zh) 心跳检测方法、装置、存储介质及终端
CN113616217A (zh) 基线漂移曲线的生成方法和装置
Hassan et al. Improved hybrid approach for enhancing protein-coding regions identification in DNA sequences
WO2019014388A1 (en) SYSTEMS, METHODS, AND MEDIA FOR INTEGRATED, EFFICIENT, REAL-TIME PROCESSING OF PHYSIOLOGICAL SIGNALS USING S-PROCESSES
Serasinghe et al. A novel low-complexity VLSI architecture for an EEG feature extraction platform
TW201236641A (en) Touch-sensitive display apparatus with pulse measuring function and method thereof
CN103829975B (zh) 超声弹性成像位移值的计算方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 13892870

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 30/06/2016)

122 Ep: pct application non-entry in european phase

Ref document number: 13892870

Country of ref document: EP

Kind code of ref document: A1