Wavenumber domain three-dimensional ultrasound full matrix imaging method based on CUDA

Abstract

The present invention discloses a wavenumber domain three-dimensional ultrasound full-matrix imaging method based on CUDA. The present invention represents the original full-matrix data in the form of a five-dimensional matrix, converts the obtained original full-matrix data into the wavenumber domain through a five-dimensional fast Fourier transform, realizes the wavenumber conversion between the original data coordinate system and the imaging coordinate system by means of Weyl's Identity, obtains the mapping of the image coordinate system through eight linear interpolation, and then obtains the three-dimensional image through a three-dimensional fast inverse Fourier transform. The present invention adopts the VTK module and uses the MIP (maximum intensity projection) ray projection method to obtain the rendered two-dimensional image, and can change the observation angle and distance of the three-dimensional image in real time. The present invention has the real-time high-efficiency, high-quality three-dimensional imaging and two-dimensional rendering capabilities for variable regions of interest, can be flexibly transplanted to computer equipment with a CUDA computing platform, and can cope with the rapid processing of large amounts of three-dimensional full-matrix data.

Description

Technical Field

[0001] The invention belongs to the technical field of industrial ultrasonic nondestructive testing, and relates to a method for realizing wavenumber domain three-dimensional ultrasonic full-matrix imaging based on CUDA. Background Art

[0002] As an important nondestructive testing technology, ultrasonic testing has shifted from nondestructive flaw detection and qualitative analysis of nondestructive testing to quantitative nondestructive evaluation, which also puts forward higher requirements on the quantitative detection characteristics of ultrasonic testing technology. With the complexity of modern industrial structure and the extreme operating conditions of industrial equipment, the requirements for the quality of workpieces and equipment are getting higher and higher. At the same time, ultrasonic testing technology is also required to be able to discover and detect more subtle defects, and there are higher requirements for the location, shape and type of detected defects. In this case, traditional conventional ultrasound can no longer meet the requirements. Phased array transducers are composed of multiple groups of small chip array elements arranged in a certain geometric manner, usually with linear arrays, planar arrays and circular arrays. With its flexible synthesis and adaptive control capabilities of the beam, phased array ultrasonic testing technology has obvious technical advantages in detection efficiency, detection range, sound beam accessibility and sensitivity, which has also enabled it to develop rapidly and be widely used in the past two decades.

[0003] In phased array ultrasonic nondestructive testing, full matrix data acquisition uses the phased array transducer to stimulate a single transducer element to emit ultrasonic waves in turn, and then uses all elements to receive and store echo signals, making full use of the signal data transmitted and received by all elements in the phased array ultrasonic transducer, including the sound wave information of the transducer element incident on the workpiece under test from different angles, with more comprehensive internal information of the workpiece under test, and a more comprehensive and accurate description of the defect location, size and shape. Especially when the defect location is unknown, phased array ultrasonic full matrix imaging does not need to repeatedly adjust the sound beam deflection angle and focusing position, and can obtain comprehensive information of unknown defects according to the original data acquisition method. This feature greatly improves the detection accuracy of phased array ultrasonic imaging. In addition, the full matrix imaging technology does not need to pre-calculate the delay focusing law, nor does it need to adjust the beam synthesis scheme according to different defect types and defect locations. It only needs to detect and collect signals for all workpieces under test according to the full matrix acquisition method, and then perform full-area focusing imaging through the host computer software, so that all defects in a certain area can be detected at one time.

[0004] Compared with traditional two-dimensional ultrasonic imaging, three-dimensional imaging has gained more and more attention due to its more comprehensive and accurate characterization of spatial defects and wider imaging area. Two-dimensional images can often only obtain information about a single cross section of a specimen, and cannot provide an overall assessment of the specimen's defects. However, in the detection of oil pipelines, assembly line inspections, and threaded fasteners, three-dimensional imaging has gradually replaced two-dimensional imaging and become an indispensable non-destructive testing method. Three-dimensional images can effectively locate the spatial position of defects in the specimen, and can render the shape and characteristics of the defects, facilitating quantitative non-destructive assessment of the damaged parts and degree of damage of the specimen.

[0005] Full-matrix 3D imaging has many conveniences, but its larger data volume and computational consumption put higher demands on post-processing capabilities, seriously affecting the efficiency and accuracy of imaging. Compared with traditional time-domain algorithms, wavenumber domain algorithms greatly improve the computational efficiency of full-matrix 3D imaging with the help of the lower computational complexity of fast Fourier transform, and the use of more accurate octal linear interpolation can effectively compensate for the phase error caused by interpolation and improve the accuracy of defect imaging in the region of interest. In order to further improve the computational efficiency on the basis of ensuring imaging accuracy, meet the requirements of real-time online detection, and make full use of the parallel computing and flexible transplantation capabilities of fast Fourier transform, the wavenumber domain algorithm can further use the GPU with CUDA platform to accelerate the calculation process, and use the VTK module to realize GPU 2D rendering, so as to achieve a complete process from raw data to image display and keyboard and mouse interaction. Summary of the invention

[0006] The purpose of the present invention is to provide a CUDA-based wavenumber domain three-dimensional ultrasound full matrix imaging method. The method of the present invention can improve the efficiency and accuracy of three-dimensional full matrix imaging, render two-dimensional maximum density projection, and achieve the goal of real-time online three-dimensional ultrasound non-destructive testing. The basic parameters of the five-dimensional fast Fourier transform are set according to the size and range of the three-dimensional image, the wavenumber domain truncation filter is set according to the center frequency and bandwidth of the array ultrasound transducer, the wavenumber conversion relationship between the original data coordinate system and the image coordinate system is obtained based on Weyl's Identity, the wavenumber domain original data is mapped to the wavenumber domain image through eight linear interpolation, and the required three-dimensional image is obtained by three-dimensional fast inverse Fourier transform. The maximum density projection of the three-dimensional image is rendered by the VTK module and displayed through an interactive window. The observation angle and distance of the three-dimensional image can be changed in real time, and the data can be refreshed in real time, which is conducive to real-time three-dimensional ultrasound non-destructive quantitative evaluation.

[0007] The technical solution adopted by the present invention is:

[0008] The present invention first provides a CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method, which comprises the following steps:

[0009] Step 1: Set up a planar array ultrasonic transducer, couple it with the isotropic homogeneous medium specimen in direct contact, use the full matrix mode to collect original three-dimensional full matrix data, determine the sound velocity of the specimen, the parameters of the planar array ultrasonic transducer, and the size range of the three-dimensional image, take the center of the array element in the lower left corner of the planar array ultrasonic transducer as the coordinate origin, take the two sides of the planar array ultrasonic transducer as the X-axis and Y-axis, and take the sound wave propagation direction perpendicular to the transducer plane as the Z-axis to establish the original data coordinate system;

[0010] Step 2: Determine the calculation length of each dimension of the five-dimensional Fourier transform, write a CUDA kernel function based on Weyl's Identity, and calculate the wave number conversion relationship between the original data coordinate system and the image coordinate system;

[0011] Step 3: Copy the original three-dimensional full matrix data in the CPU memory to the GPU video memory, rearrange the three-dimensional full matrix data into five-dimensional full matrix data in the CUDA kernel function, and the data index is composed by dimension; use the row-column algorithm to decompose the five-dimensional Fourier transform, change the index of the data in the storage space according to the coordinate relationship of the matrix transposition, and convert the five-dimensional Fourier transform into calculating multiple two-dimensional and three-dimensional Fourier transforms; use the fast five-dimensional Fourier transform to convert the original data into the wave number domain, and perform rectangular window truncation on the frequency dimension;

[0012] Step 4: For each set of three-dimensional data in the emission dimension of the five-dimensional matrix in the wavenumber domain, the wavenumber conversion relationship between the original data coordinate system and the image coordinate system in step 2 is used as the interpolation index, and several adjacent points in the spatial dimension are selected as references. The Stolt mapping of the five-dimensional matrix in the wavenumber domain under the three-dimensional image is calculated using the interpolation algorithm;

[0013] Step 5: For the wavenumber domain Stolt mapping calculated in step 4, extract each set of three-dimensional data, align and superimpose them respectively, obtain a three-dimensional wavenumber domain image matrix, perform a three-dimensional fast inverse Fourier transform, transform the image matrix from the wavenumber domain to the spatial domain, generate a three-dimensional ultrasound full matrix image, and truncate the part outside the region of interest;

[0014] Step 6: For the three-dimensional ultrasound full matrix image generated in step 5, copy it from the GPU video memory to the CPU memory, push it into the data processing pipeline of the VTK module, use maximum density projection to render it into a two-dimensional image and display it in the window. It can also update the data in real time and interactively change the viewing position and focus distance.

[0015] As a preferred embodiment of the present invention, the test piece described in step 1 is a uniform isotropic medium, the sound velocity of which is set to c, and the number of array elements of the surface array ultrasonic transducer used is distributed as follows: The discrete subscripts w1 and w2 represent the number of array elements of the array ultrasonic transducer in the X-axis and Y-axis directions, respectively, and the center pitch of the array elements in the X-axis and Y-axis directions are The specimen and the imaging area are placed on the positive half axis of the Z axis, and the number of pixels in the imaging area in the X, Y and Z axis directions is set to N respectively. x ,N y ,N z , the pixel spacing is d x ,d y ,d z ; The origin of the imaging area coincides with the origin of the original data coordinate system; the original three-dimensional full matrix ultrasound data obtained is recorded as e raw (t0,v,u).

[0016] As a preferred embodiment of the present invention, the step 2 is specifically as follows:

[0017] The required point values ​​for each dimension of the five-dimensional Fourier transform are obtained from the depth, width, and pixel spacing, which are N t , The discrete subscript t represents time, v1 and v2 represent the components of the receiving array element in the X and Y axis directions, and u1 and u2 represent the components of the transmitting array element in the X and Y axis directions. When full matrix acquisition is used for imaging,

[0018] Through Weyl's Identity and the forward propagation model including Green's function, the specific expression of Weyl's Identity in three-dimensional space is:

[0019]

[0020] Among them, G(ω,x,y,z) is the Green function of plane wave decomposition in three-dimensional space, k is the spatial wave number, k x ,k y are the components of the spatial wave number k on the X and Y axes, ω is the angular frequency, and x, y, z are the spatial displacements in the X, Y and Z axes respectively; according to the simplified forward propagation model of the imaging system, the transmit-receive signal pair, that is, the received data, is e(t,v1,v2,u1,u2),

[0021] Its spectrum E(ω,v1,v2,u1,u2) with respect to time can be calculated using the following relationship:

[0022]

[0023] Among them, function G is the Green function, f(x, y, z) is the defect scattering point distribution in the forward propagation model, that is, the three-dimensional image to be obtained. Substitute Weyl's Identity into the forward propagation model formula and rewrite the relationship with v1, v2, u1, u2 as variables:

[0024]

[0025]

[0026] Among them, k v1 , are the wave number components of the receiving array element in the X and Y axis directions and the wave number components of the transmitting array element in the X and Y axis directions, respectively. x ,k y ,k z ) is the wave number domain form of the three-dimensional image f(x,y,z). Therefore, the five-dimensional Fourier transform form of e(t,v1,v2,u1,u2) It can be expressed as the following relationship:

[0027]

[0028] The following substitution relationship is included:

[0029]

[0030]

[0031]

[0032] Rewrite the above formula as k x ,k y ,k z To represent k, In order to distinguish it from the spatial wave number of the original data, k is used here. g ,k v1g ,k v2g replace:

[0033] k v1g =k x -k u1g ,

[0034] k v2g =k y -k v1g ,

[0035]

[0036] where k x ,k y ,kz are the components of the spatial wave number k on the X, Y, and Z axes, respectively, where the discrete subscripts x, y, and z represent the X, Y, and Z axes, respectively, and k v1 , They are the wave number components of the receiving array element in the X and Y axis directions and the wave number components of the transmitting array element in the X and Y axis directions respectively.

[0037] For given array probe parameters and initial imaging area, its spatial wave number and wave number components in each direction are known, as expressed by the following relationship:

[0038]

[0039]

[0040]

[0041]

[0042]

[0043]

[0044]

[0045]

[0046]

[0047] f s is the sampling frequency of receiving the full matrix signal. After Fourier transform according to step 3, rectangular truncation is performed according to the probe bandwidth, and the minimum and maximum frequencies of the truncation are set to f min ,f max , the corresponding angular frequency range ω∈[2πf min ,2πf max ]. Therefore, the angular frequency after truncation is:

[0048]

[0049] According to step 3 and step 4, in order to make full use of the original data of the interpolation mapping, the spatial wave number k of the three-dimensional image in the Z-axis direction z Assume the following relationship:

[0050]

[0051] Write a CUDA kernel function to calculate the wave number in the original data coordinate system k, and image coordinate system k x ,k y ,k zThe conversion relationship is used as a reference index for wavenumber domain imaging calculation.

[0052] As a preferred embodiment of the present invention, the step three is specifically as follows: the three-dimensional original data e raw (t0,v,u) is converted into five-dimensional data e(t,v1,v2,u1,u2), where the index t0 represents the time in the original data, and the index v,u represent the receiving and transmitting array elements of the full matrix acquisition respectively. It is expressed in discrete matrix form as e raw The data index relationship between [t0][v][u] and e[t][v1][v2][u1][u2] is as follows:

[0053] t=t0

[0054]

[0055]

[0056]

[0057]

[0058] Among them, the function floor represents the function of rounding down; the five-dimensional Fourier transform is decomposed using the row-column algorithm, converted into multiple three-dimensional and two-dimensional Fourier transforms respectively, and the Fourier transform results and intermediate steps are recorded as E n (t,v1,v2,u1,u2) The calculation process is as follows:

[0059] (4) Secondary Dimension Point Fast Fourier Transform:

[0060]

[0061] (5) Transpose the data matrix

[0062]

[0063] (6) Sub-3D N t , Point Fast Fourier Transform:

[0064]

[0065] Truncating the data, retaining ω∈[2πf min ,2πf max ], let its points be

[0066] As a preferred embodiment of the present invention, the interpolation algorithm described in step 4, namely the Stolt inverse transform, needs to lock The forward propagation model formula in step 2 can be rewritten as:

[0067]

[0068] Since defect imaging detection focuses more on relative amplitude, the proportionality coefficient is ignored. Substitute the wave number parameters pre-calculated in step 2. Since the spatial wave number is nonlinear in the conversion between the original data and the image data, and lock k u1 ,k u2 After that, only ω,k v1 ,k v2 For the three parameters, the present invention adopts octal linear interpolation as the interpolation algorithm and completes the calculation through the corresponding CUDA kernel function.

[0069] As a preferred embodiment of the present invention, in step 5, the wave number domain is superimposed along all the transmitting array elements, and a three-dimensional fast inverse Fourier transform is performed, and its expression is as follows:

[0070]

[0071]

[0072] Among them, f(z,x,y) is the calculated three-dimensional ultrasound image, and the part outside the region of interest is truncated to obtain the image I(z,x,y) within the region of interest.

[0073] As a preferred embodiment of the present invention, the data processing pipeline pushed into the VTK module described in step 6 is rendered into a two-dimensional image using maximum density projection, specifically comprising the following steps:

[0074] a: Initialize VTK module parameters;

[0075] b: Filtering operation on 3D image data;

[0076] c: Rendering the 3D image data into a 2D image using maximum intensity projection;

[0077] d: Set up the VTK interactive window to display and refresh the two-dimensional image in real time.

[0078] The beneficial effects of the present invention are mainly manifested in:

[0079] (1) The present invention utilizes fast Fourier transform to convert conventional time-domain delay superposition into wavenumber domain processing, which greatly reduces the computational complexity of previous time-domain algorithms, thereby improving the efficiency and accuracy of three-dimensional ultrasound imaging and achieving real-time high frame rate three-dimensional ultrasound full-matrix imaging.

[0080] (2) The present invention adopts an efficient and concise wavenumber domain calculation formula, which can calculate the values ​​of all pixel points in the image as a whole. The formula is highly stable and the method is simple and easy to understand. It avoids the inefficiency of conventional time domain algorithms that delay and superimpose calculations of each pixel point separately, thereby improving the speed of three-dimensional ultrasound full-matrix imaging in principle.

[0081] (3) The present invention adopts the spectrum processing technology in the wave number domain, which reduces the calculation time of means such as filtering to improve image quality, thereby effectively reducing the calculation time and obtaining a higher imaging effect at a lower calculation cost.

[0082] (4) The present invention adopts an interpolation method that is independently written and modified to solve the problem that it is difficult to implement complex interpolation methods in software, thereby effectively solving the software implementation problem of three-dimensional wavenumber domain ultrasonic imaging, and its calculation speed is relatively fast, which also reduces the calculation consumption.

[0083] (5) The present invention utilizes the GPU parallel processing capability of CUDA to further improve the efficiency of three-dimensional imaging. Since the proposed fast Fourier transform conversion into wavenumber domain processing, spectrum processing technology and interpolation method are all implemented in the present invention using the CUDA platform method, the main calculation is performed in the GPU, which solves the low efficiency of serial calculation in conventional CPU methods and greatly saves imaging time.

[0084] (6) The algorithm of the present invention is extremely convenient and portable. It uses a computing device with a CUDA platform, which avoids the problem of relatively long code burning and code update cycles compared to ordinary FPGAs (field programmable gate arrays), facilitates error correction and personal customization, and enables three-dimensional real-time ultrasonic imaging to be achieved in complex and flexible on-site environments.

[0085] (7) The algorithm of the present invention not only provides three-dimensional ultrasound images, but also can further render two-dimensional mapping images for display on the screen after obtaining the three-dimensional ultrasound images. It also has the ability to update the image source and camera viewing angle in real time, and the rendering window can interact in real time. This avoids the problem that conventional methods only calculate three-dimensional images and cannot provide users with intuitive and real-time imaging detection capabilities, which greatly benefits users' ability to read and distinguish images on site and in post-processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1 It is an overall flow chart of the method of the present invention.

[0087] Figure 2 It is a schematic diagram of the installation position of the 3D detection model.

[0088] Figure 3 It uses the VTK module to render and display a two-dimensional image flow chart in real time.

[0089] Figure 4It uses three-dimensional real-time ultrasound imaging and real-time rendering of two-dimensional image results. DETAILED DESCRIPTION

[0090] The present invention will be further described below in conjunction with the accompanying drawings.

[0091] Embodiments of the present invention are as follows:

[0092] like Figure 1 As shown, the operation flow of the present invention using a CUDA-based wavenumber domain three-dimensional ultrasound full matrix imaging method can be divided into the following steps:

[0093] Step 1: Set up the area array phased array probe, couple it with the isotropic homogeneous medium specimen in direct contact, and use the full matrix mode to collect raw data. Figure 2 As shown, the ultrasonic transducer is an 8×8 square array probe, and the center spacing of each array element in the X and Y axis directions is d x =1.5mm,d y =1.5mm, each array element is square, and the gap between the array elements is s=0.1mm in the X-axis and Y-axis directions respectively. The specimen is selected as metal aluminum, whose sound velocity is c=6350m / s. The three-dimensional imaging area is below the array probe, and its origin coincides with the origin of the array probe, and the depth is h z The widths in the X and Y directions are h x and h y , the pixel spacing in each direction is δ x ,δ y ,δ z , thus calculating the number of pixels N in each direction x =ceil(h x / δ x ),N y =ceil(h y / δ y ),N z =ceil(h z / δ z ).

[0094] Step 2: Calculate the required point values ​​for each dimension of the five-dimensional Fourier transform from the depth, width, and pixel spacing, which are N t , Through Weyl's Identity and the forward propagation model including Green's function, the coordinate axis wavenumber conversion relationship between the original coordinate system and the image coordinate system in the wavenumber domain is derived:

[0095]

[0096]

[0097]

[0098] The specific implementation requires changing the composition of the variables and rewriting the above formula into k x ,k y ,k z To represent k, In order to distinguish it from the spatial wave number of the original data, k is used here. g ,k v1g ,k v2g Instead, the subscripts g, v1g, and v2g represent the spatial wavenumbers in the imaging coordinate system, the wavenumbers in the X-axis direction, and the wavenumbers in the Y-axis direction:

[0099] k v1g =k x -k u1g ,

[0100] k v2g =k y -k v1g ,

[0101]

[0102] Due to the band-limited characteristics of the ultrasonic signal, a rectangular cutoff is made according to the probe bandwidth, and the minimum and maximum cutoff frequencies are set to f min ,f max , the corresponding angular frequency range ω∈[2πf min ,2πf max ]. Therefore, the angular frequency after truncation is:

[0103]

[0104] Step 3: Transform the original full-matrix three-dimensional data into a five-dimensional matrix, and convert the data to the wavenumber domain by performing multiple two-dimensional and three-dimensional Fourier transforms instead of direct five-dimensional Fourier transforms.

[0105] The three-dimensional original data raw (t0,v,u) is converted to five-dimensional data e(t,v1,v2,u1,u2), expressed as a discrete matrix form e raw The data index relationship between [t0][v][u] and e[t][v1][v2][u1][u2] is as follows:

[0106] t=t0

[0107]

[0108]

[0109]

[0110]

[0111] The five-dimensional Fourier transform is decomposed using the row-column algorithm and converted into multiple three-dimensional and two-dimensional Fourier transforms respectively. The Fourier transform result and the intermediate steps are recorded as E n (t,v1,v2,u1,u2) The calculation process is as follows:

[0112] (7) Secondary Dimension Point Fast Fourier Transform:

[0113]

[0114] (8) Transpose the data matrix

[0115]

[0116] (9) Sub-3D N t , Point Fast Fourier Transform:

[0117]

[0118] Truncating the data, retaining ω∈[2πf min ,2πf max ], let its points be

[0119] Step 4: The calculated wavenumber domain data is obtained through k g ,k v1g ,k v2g For the three-dimensional data matrix corresponding to each emission, eight linear interpolation calculations are performed respectively and converted to the corresponding image coordinate system. The interpolation index relationship in the X, Y and Z axis directions is as follows:

[0120] (1) All data points are read in parallel. First, the point with the smallest coordinates on the three directions of X, Y, and Z in the cube required for interpolation is obtained:

[0121]

[0122]

[0123] ind z =floor((k z -min(k)) / (k(1)-k(0)))

[0124] (2) Calculate the remaining 7 points of the cube from the reference point and truncate the points that exceed The coordinate index of is the boundary value.

[0125] (3) According to the locked Substitute them separately to calculate the common product coefficient bord in each coordinate axis direction used for interpolation a ,bord b ,bord c ,a,b,c:

[0126] bord a =(1-ind x ) / N x

[0127] bord b =(1-ind y ) / N y

[0128]

[0129] a=k v1g -ind x

[0130] b=k v2g -ind y

[0131] c=k g -ind z

[0132] (4) Lock Calculate the eight linear interpolation results from the common coefficients and wavenumber domain data:

[0133]

[0134] Step 5: Align and superimpose the wavenumber domain data in the image coordinate system after interpolation corresponding to each group of transmissions obtained in step 4, and the sum is a set of three-dimensional data matrices. The matrix is ​​subjected to a three-dimensional fast inverse Fourier transform to obtain a three-dimensional ultrasound full matrix image. According to the size of the region of interest N in each direction x_roi ,N y_roi ,N z_roi Truncates the image.

[0135] In step 5, the matrix undergoes a three-dimensional fast inverse Fourier transform, and its expression is as follows:

[0136]

[0137]

[0138] Among them, f(z,x,y) is the calculated three-dimensional ultrasound image, and the part outside the region of interest is truncated to obtain I(z,x,y).

[0139] Step 6: Copy the data from the GPU video memory to the internal memory, use the maximum density projection of the VTK module to render the two-dimensional image, reduce the grid spacing during the data mapping process, improve the image quality, set the camera view and image window, and increase the window interaction function. In order to improve the rendering effect during the rendering process, use a 10x upsampling factor. The process of using the VTK module to render and display a two-dimensional image in real time in this embodiment is as follows: Figure 3 As shown, specifically including:

[0140] a: Initialize module parameters, set VTK data origin and grid spacing, set VTK volume mapping memory space and sampling spacing, set VTK volume properties, including color and opacity, etc.

[0141] b: Intermediate processing of 3D image data, performing upsampling, downsampling and filtering operations to obtain the processed 3D image.

[0142] c: The processed 3D image data is rendered into 2D image data using maximum density projection. The volume mapping and attributes need to be added to the VTK volume rendering variable. The VTK rendering variable and camera perspective are set, and the 3D image data is imported into the rendering variable to obtain the 2D image data.

[0143] d: Set the VTK interactive window variables, import the two-dimensional image data, and display and refresh the window image in real time according to the changes of the original full matrix data. The camera perspective can also be changed in real time.

[0144] Design a square metal specimen with five hole-type defects arranged in an oblique direction. The five holes are of the same size and are arranged equidistantly in the oblique direction. t =1024, For example, the number of pixels in the three-dimensional imaging area is selected as N x =64,N y =64,N z =512, after the five-dimensional Fourier transform, the minimum and maximum frequencies of the frequency dimension are truncated to f min =3MHz,f max =7MHz. The approximate range where the defect is located is selected as the region of interest, and its size is N x_roi =64,N y_roi =64,N z_roi =300, through wave number domain calculation and VTK rendering, the imaging results are as follows Figure 4As shown, the intensity of the pixels can be used completely and hierarchically to show the defects of the five designed holes in the image, and the square bottom surface of the specimen can also be completely calculated. The image quality is high and there is no sense of stripes.

[0145] Due to the low algorithm complexity of fast Fourier transform and the parallel computing capability of GPU, large-scale real-time high-frame-rate three-dimensional ultrasound full-matrix imaging, rendering and display can be achieved.

[0146] The above-mentioned embodiments only express several implementation methods of the present invention, and the description thereof is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. It should be pointed out that, for a person of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the attached claims.

Claims

1. A CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method, characterized in that: The steps include: Step 1: Set up a planar array ultrasonic transducer, couple it with the isotropic homogeneous medium specimen in direct contact, use the full matrix mode to collect original three-dimensional full matrix data, determine the sound velocity of the specimen, the parameters of the planar array ultrasonic transducer, and the size range of the three-dimensional image, take the center of the array element in the lower left corner of the planar array ultrasonic transducer as the coordinate origin, take the two sides of the planar array ultrasonic transducer as the X-axis and Y-axis, and take the sound wave propagation direction perpendicular to the transducer plane as the Z-axis to establish the original data coordinate system; Step 2: Determine the calculation length of each dimension of the five-dimensional Fourier transform, write a CUDA kernel function based on Weyl's Identity, and calculate the wave number conversion relationship between the original data coordinate system and the image coordinate system; The specific step 2 is: calculate the required point values ​​of each dimension of the five-dimensional Fourier transform from the depth, width and pixel spacing, which are respectively The discrete subscript t represents time, v1, v2 represent the components of the receiving array element in the X and Y axis directions, and u1, u2 represent the components of the transmitting array element in the X and Y axis directions; When imaging using full matrix acquisition, Through Weyl's Identity and the forward propagation model including Green's function, the specific expression of Weyl's Identity in three-dimensional space is: Among them, G(ω,x,y,z) is the Green function of plane wave decomposition in three-dimensional space, k is the spatial wave number, k x ,k y are the components of the spatial wave number k on the X and Y axes, ω is the angular frequency, and x, y, z are the spatial displacements in the X, Y and Z axes respectively; according to the simplified forward propagation model of the imaging system, the transmit-receive signal pair, that is, the received data, is e(t,v1,v2,u1,u2), Its spectrum E(ω,v1,v2,u1,u2) with respect to time is calculated using the following relationship: Among them, function G is the Green function, f(x, y, z) is the defect scattering point distribution in the forward propagation model, that is, the three-dimensional image to be obtained; Substitute Weyl's Identity into the forward propagation model formula and rewrite the relationship with v1, v2, u1, u2 as variables: in, are the wave number components of the receiving array element in the X and Y axis directions and the wave number components of the transmitting array element in the X and Y axis directions, respectively. x ,k y ,k z ) is the wave number domain form of the three-dimensional image f(x,y,z); therefore, the five-dimensional Fourier transform form of e(t,v1,v2,u1,u2) It is expressed as the following relationship: The following substitution relationship is included: Rewrite the above formula as To express In order to distinguish it from the spatial wave number of the original data, k is used here. g ,k v1g ,k v2g replace: where k x ,k y ,k z are the components of the spatial wave number k on the X, Y, and Z axes, respectively, where the discrete subscripts x, y, and z represent the X, Y, and Z axes, respectively. are the wave number components of the receiving array element in the X and Y axis directions and the wave number components of the transmitting array element in the X and Y axis directions respectively; Step 3: Copy the original three-dimensional full matrix data in the CPU memory to the GPU video memory, rearrange the three-dimensional full matrix data into five-dimensional full matrix data in the CUDA kernel function, and the data index is composed by dimension; use the row-column algorithm to decompose the five-dimensional Fourier transform, change the index of the data in the storage space according to the coordinate relationship of the matrix transposition, and convert the five-dimensional Fourier transform into calculating multiple two-dimensional and three-dimensional Fourier transforms; use the fast five-dimensional Fourier transform to convert the original data into the wave number domain, and perform rectangular window truncation on the frequency dimension; Step 4: For each set of three-dimensional data in the emission dimension of the five-dimensional matrix in the wavenumber domain, the wavenumber conversion relationship between the original data coordinate system and the image coordinate system in step 2 is used as the interpolation index, and several adjacent points in the spatial dimension are selected as references. The Stolt mapping of the five-dimensional matrix in the wavenumber domain under the three-dimensional image is calculated using the interpolation algorithm; Step 5: For the wavenumber domain Stolt mapping calculated in step 4, extract each set of three-dimensional data, align and superimpose them respectively, obtain a three-dimensional wavenumber domain image matrix, perform a three-dimensional fast inverse Fourier transform, transform the image matrix from the wavenumber domain to the spatial domain, generate a three-dimensional ultrasound full matrix image, and truncate the part outside the region of interest; Step 6: For the three-dimensional ultrasound full matrix image generated in step 5, copy it from the GPU video memory to the CPU memory, push it into the data processing pipeline of the VTK module, use maximum density projection to render it into a two-dimensional image and display it in the window. It can also update the data in real time and interactively change the viewing position and focus distance.

2. The CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method according to claim 1, characterized in that: The test piece described in step 1 is a uniform isotropic medium, the sound velocity of which is set to c, and the number of array elements of the array array ultrasonic transducer used is N. w1 ×N w2 , where the discrete subscripts w1 and w2 represent the number of array elements of the array ultrasonic transducer in the X-axis and Y-axis directions, respectively, and the pitch of the array element centers in the X-axis and Y-axis directions are d w1 ,d w2 The specimen and the imaging area are placed on the positive half axis of the Z axis, and the number of pixels in the imaging area in the X, Y and Z axis directions is set to N respectively. x ,N y ,N z , the pixel spacing is d x ,d y ,d z ; The origin of the imaging area coincides with the origin of the original data coordinate system; the original three-dimensional full matrix ultrasound data obtained is recorded as e raw (t0,v,u).

3. The CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method according to claim 1, characterized in that: The step three is specifically as follows: the three-dimensional original data e raw (t0,v,u) is converted into five-dimensional data e(t,v1,v2,u1,u2), where the index t0 represents the time in the original data, and the index v,u represent the receiving and transmitting array elements of the full matrix acquisition, respectively, expressed in discrete matrix form as e raw The data index relationship between [t0][v][u] and e[t][v1][v2][u1][u2] is as follows: t=t0 Among them, the function floor represents the function of rounding down; the five-dimensional Fourier transform is decomposed using the row-column algorithm, converted into multiple three-dimensional and two-dimensional Fourier transforms respectively, and the Fourier transform results and intermediate steps are recorded as E n (t,v1,v2,u1,u2) The calculation process is as follows: (1) Secondary Dimension Point Fast Fourier Transform: (2) Transpose the data matrix (3) Sub-3D Point Fast Fourier Transform: Truncating the data, retaining ω∈[2πf min ,2πf max ], let its points be 4. The CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method according to claim 1, characterized in that: The interpolation algorithm in step 4 is octal linear interpolation.

5. The CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method according to claim 1, characterized in that: In the step 5, all the transmitting array elements are superimposed in the wave number domain and a three-dimensional fast inverse Fourier transform is performed, and the expression is as follows: Among them, f(z,x,y) is the calculated three-dimensional ultrasound image, and the part outside the region of interest is truncated to obtain the image I(z,x,y) within the region of interest.

6. The CUDA-based wavenumber domain three-dimensional ultrasound full-matrix imaging method according to claim 1, characterized in that: The data processing pipeline pushed into the VTK module described in step 6 is rendered into a two-dimensional image using maximum density projection, and specifically includes the following steps: a: Initialize VTK module parameters; b: Filtering operation on 3D image data; c: Rendering the 3D image data into a 2D image using maximum intensity projection; d: Set up the VTK interactive window to display and refresh the two-dimensional image in real time.

Images