CN111248858A - Photoacoustic tomography reconstruction method based on frequency domain wavenumber domain - Google Patents

Abstract

The invention discloses a photoacoustic tomography reconstruction method based on a frequency domain wavenumber domain, which comprises the following steps of: the method comprises the steps of obtaining radio frequency data and system parameters, performing time gain compensation, performing three-dimensional Fourier transform, performing interpolation processing, performing three-dimensional inverse Fourier transform, demodulating envelope, adjusting contrast, outputting an optical acoustic image, and optimizing the optical acoustic image. The invention carries out three-dimensional reconstruction of the photoacoustic signal in the frequency domain and the wave number domain, reduces the influence of incomplete signal acquisition under a limited angle, and simultaneously considers the influence of sound velocity nonuniformity. The reconstruction method provided by the invention has the advantages of high resolution of the reconstructed image, high signal-to-noise ratio, less artifacts and high reconstruction calculation speed.

Description

A Reconstruction Method of Photoacoustic Tomography Based on Frequency Domain and Wavenumber Domain

technical field

The invention belongs to an imaging reconstruction method, in particular to a photoacoustic tomographic imaging reconstruction method.

Background technique

Photoacoustic imaging is an emerging imaging technique that offers potential perspectives and opportunities for clinical imaging. Photoacoustic tomography is an important branch of photoacoustic imaging, and its advantage lies in the deep imaging depth, which can reach several centimeters, which meets the basic needs of clinical imaging. At present, photoacoustic tomography based on annular array acoustic transducer (full angle detection) has been used for breast examination. However, the annular acoustic detector is limited by its inherent shape and structure, and it is difficult to be used for photoacoustic imaging and imaging of more parts. clinical examination. The acoustic transducer based on hand-held (limited angle detection) has the advantage of being flexible and portable, and at the same time, it is easy to combine with ultrasound imaging commonly used in clinical practice, which is the development trend of photoacoustic tomography. The most notable feature of photoacoustic tomography that is different from other photoacoustic imaging techniques is that photoacoustic tomography needs to reconstruct the collected radio frequency data to obtain photoacoustic images. Currently, the commonly used reconstruction methods include time inversion method and filtered back projection method. , these two methods have good effect on photoacoustic tomography imaging with full-angle detection; however, for hand-held acoustic transducers, because only signals from a limited angle are collected, the above two methods are difficult to achieve accurate reconstruction imaging. On the other hand, the two-dimensional matrix acoustic transducer is an emerging handheld acoustic transducer that can directly perform three-dimensional imaging, and there is no method that can reconstruct the two-dimensional matrix acoustic transducer. Therefore, a new type of photoacoustic reconstruction algorithm is urgently needed to realize accurate reconstruction imaging under the limited angle detection of photoacoustics based on two-dimensional matrix acoustic transducers.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a method for reconstructing photoacoustic tomography based on the frequency domain and wavenumber domain, which is suitable for a two-dimensional matrix acoustic transducer, and can perform optical imaging in the frequency domain and the wavenumber domain. The three-dimensional reconstruction of the acoustic signal also reduces the influence of incomplete signal acquisition under a limited angle, and the reconstruction calculation amount is small. In addition, because it is not based on the spatial projection of the time domain signal, the influence of projection artifacts is greatly reduced.

In order to achieve the above purpose, the photoacoustic tomography reconstruction method based on the frequency domain and wavenumber domain of the present invention includes the following steps:

(1) Acquire the radio frequency data R (x, y, z=0, t) collected by the two-dimensional matrix acoustic transducer, as well as the system parameters of the two-dimensional matrix acoustic transducer and the sound velocity distribution V ( z'). The two-dimensional matrix acoustic transducer is an ultrasonic transducer based on a linear array, the radio frequency data R(x, y, z=0, t) is collected in real time or saved in advance, and x, y are along the two-dimensional matrix The position coordinates of the acoustic transducer, and x, y are perpendicular to each other, z is the position coordinate of the vertical two-dimensional matrix acoustic transducer direction, t is the coordinate representing the time sequence, and the system parameters of the two-dimensional matrix acoustic transducer is the sampling frequency fs of the two-dimensional matrix acoustic transducer, the number n of the array elements of the two-dimensional matrix acoustic transducer, and the array element spacing pt of the two-dimensional matrix acoustic transducer.

(2) Perform time gain compensation on the collected radio frequency data R (x, y, z=0, t) to obtain compensated radio frequency data s (x, y, z=0, t). The specific method of time gain compensation is to multiply the value on each column of the radio frequency data R(x, y, z=0, t) by the weighting function Y(t), where Y(t)=at+1. That is, s(x, y, z=0, t)=R(x, y, z=0, t)×(at+1), 1<a<5; since a is greater than 1, it increases with time , the RF data will be processed accordingly to counteract the loss of the acoustic signal during propagation.

因此对Φ(k_x，k_y，z＝0，f)进行插值即可得到Φ′(k_x，k_y，k_z，t＝0)；此处需要强调的是频域数据Φ(k_x，k_y，z＝0，f)与原始采集到的射频数据互为傅里叶变换对，而空间波数域数据Φ′(k_x，k_y，k_z，t＝0)与待重建的图像互为傅里叶变换对。作为优选的，插值方法是Φ′(k_x，k_y，k_z，t＝0)＝∫Φ(kx，ky，z＝0，f)M(k_x，k_y，k_z，f)df，其中：

(3) Perform three-dimensional Fourier transform on the compensated radio frequency data s(x, y, z=0, t) according to x, y and t to obtain frequency domain data Φ(k _x , _ky , z=0, f ), and then interpolate Φ((k _x , _ky , z=0, f) to obtain spatial wavenumber domain data Φ′(k _x , _ky , k _z , t=0), where k _x is the direction along x The wave number in the direction, ky is the wave number along the y direction, k _z is the wave number along the _z direction, and f is the frequency. It can be known from the wave equation that the frequency domain of the z=0 plane (where the two-dimensional matrix acoustic transducer is located) The data Φ(k _x , _ky , z=0, f) differs by one from the spatial wavenumber domain data Φ′(k _x , _ky , k _z , t=0) of the initial position (t=0) in the detection plane Phase factor

Therefore, Φ(k _x , _ky , z=0, f) can be interpolated to obtain Φ′(k _x , _ky , k _z , t=0); what needs to be emphasized here is the frequency domain data Φ(k _x , _ky , z=0, f) and the originally collected radio frequency data are Fourier transform pairs, while the spatial wavenumber domain data Φ′ (k _x , _ky , k _z , t=0) and the data to be reconstructed The images are Fourier transform pairs of each other. Preferably, the interpolation method is Φ'(k _x , _ky , k _z , t=0)=∫Φ(kx, ky, z=0, f)M(k _x , _ky , k _z , f) df, where:

(4), perform three-dimensional inverse Fourier transform on spatial frequency domain data Φ' (k _x , _ky , k _z , t=0) to obtain spatial amplitude data s' (x, y, z, t=0) .

(5) Demodulate the spatial amplitude data s' (x, y, z, t=0) to obtain the envelope, and then adjust the contrast, and finally obtain the photoacoustic reconstructed image s*, in which the demodulation takes the envelope The method is the quadrature demodulation method or the Hilbert demodulation method, and the contrast adjustment method can be a logarithmic compression method, or the data obtained after de-enveloping can be raised to the power of n, where 0<n<2.

(6) Optimizing the sound velocity distribution v(z'), and substituting the optimized V ₂ (z') data into steps (3) to (5) to obtain a photoacoustic reconstructed image S with the highest tissue contrast CTR ₂ *;

The specific steps for optimizing the sound velocity distribution V(z') are:

(1), set the optimization times counter t=0, set the maximum optimization times T, the value range of T is 3-10, and randomly generate N sound velocity distributions V ₁ (z′) as the initial sample set K ₁ , the value of N is The value range is 300-1000, where the range of sound speed in each sound speed distribution is randomly within 1350-1650m/s;

(2) Reconstructing the photoacoustic image for each sound velocity distribution V ₁ (z′) in the above steps, and the reconstruction steps are the steps (1) to (5) in the above-mentioned photoacoustic tomography reconstruction method. According to the reconstruction method described above, N photoacoustic reconstructed images Si* are obtained, and the CTR of each _{photoacoustic} reconstructed image _Si * is calculated. The specific calculation method of CTR is the grayscale of the target object and the background in the _{photoacoustic} reconstruction image Si*. ratio of values;

(3) Sort the N calculated _CTRs in descending order of numerical value, select the sound velocity distribution corresponding to the first 30% larger CTRs to form the sample set Ka, and select the sound velocity distribution corresponding to the first 30%-60% larger CTRs to form The number of samples in the sample set K _b , K _a and K _b are all 0.3N;

(4) Randomly select a sound speed distribution sample V a (z') from Ka, randomly select a sound speed distribution sample V _{b (} z' _{) from K b ,} and take the mean value of the two sound speed distribution samples as the fused Sound velocity distribution V _c (z');

(5) Repeat step (4) to obtain 0.2N fused sound velocity distributions V _c (z′), and compare the unselected sound velocity distributions in step (4 _{) Ka and K b with} the fused sound velocity distribution V _c (z') together constitute a sample set K ₂ , and the number of samples in K ₂ is 0.4N;

(6), repeat steps (2)-(6) with K ₂ as the initial sample set, continuously update the sample set of the sound velocity distribution;

(7) When t=T, the optimization of the sound velocity distribution V(z′) is completed, and the sample with the largest CTR in the final sample set is selected as the final sound velocity distribution V ₂ (z′).

Compared with the prior art, the present invention has the following beneficial effects:

For the photoacoustic tomography reconstruction based on the hand-held two-dimensional matrix acoustic transducer, the reconstruction method proposed in the present invention has the advantages of high reconstructed image resolution, high signal-to-noise ratio, less artifacts, and fast reconstruction calculation speed.

Description of drawings

Fig. 1 is the flow chart of the present invention;

Figure 2 is a schematic diagram of a two-dimensional matrix acoustic transducer acquisition;

Fig. 3 is the reconstruction result of point sound source RF data based on simulation;

FIG. 4 is a two-dimensional cross-section of the reconstruction result of the vascular radiofrequency data based on simulation;

FIG. 5 is a two-dimensional cross-section of the reconstruction result of the RF data of the blood vessel model based on the real detection.

Detailed ways

In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

As shown in FIG. 1 , it is a flow chart of a photoacoustic tomography reconstruction method based on the frequency domain and wavenumber domain of the present invention.

Step (1), obtain the radio frequency data R (x, y, z=0, t) collected by the two-dimensional matrix acoustic transducer, as well as the system parameters of the two-dimensional matrix acoustic transducer and the sound velocity distribution V of the object to be measured (z'). The two-dimensional matrix acoustic transducer is an ultrasonic transducer based on a linear array, and the radio frequency data R(x, y, z=0, t) is collected in real time or saved in advance; as shown in Figure 2, the two Schematic diagram of the acquisition of the 2D matrix acoustic transducer, x, y are the position coordinates along the 2D matrix acoustic transducer, and x, y are perpendicular to each other, z is the position coordinate perpendicular to the direction of the 2D matrix acoustic transducer, t is the coordinate representing the time sequence; the system parameters of the two-dimensional matrix acoustic transducer are the sampling frequency fs of the two-dimensional matrix acoustic transducer, the number of elements n of the two-dimensional matrix acoustic transducer, the two-dimensional matrix acoustic transducer The array element spacing pt and the sound velocity distribution v ² (z'). The width of the reconstructed image is pt×n, and the height of the reconstructed image is determined by the length of the radio frequency data, the sampling frequency fs of the two-dimensional matrix acoustic transducer, and the speed of sound.

In step (2), time gain compensation is performed on the collected radio frequency data R(x, y, z=0, t) to obtain compensated radio frequency data s (x, y, z=0, t). The specific method of time gain compensation is to multiply the value on each column of the radio frequency data R(x, y, z=0, t) by the weighting function Y(t), Y(t)=2t+1. That is, s(x, y, z=0, t)=R(x, y, z=0, t)×(2t+1); as time increases, the radio frequency data will be counteracted accordingly to cancel The loss of an acoustic signal during propagation.

Step (3), perform two-dimensional Fourier transform on the compensated radio frequency data s(x, y, z=0, t) according to x and t to obtain frequency domain data Φ(k _x , _ky , z=0, f), and then interpolate Φ((k _x , _ky , z=0, f) to obtain spatial wavenumber domain data Φ′(k _x , _ky , k _z , t=0), where k _x is the edge The wave number in the x direction, ky is the wave number in the y direction, k _z is the wave number in the _z direction, and f is the frequency.

因此对Φ(k_x，k_y，z＝0，f)进行插值即可得到Φ′(k_x，k_y，k_z，t＝0)。It can be known from the wave equation that the frequency domain data Φ(k _x , _ky , z=0, f) of the z=0 plane (where the two-dimensional matrix acoustic transducer is located) is related to the initial position in the detection plane (t=0 ) of the spatial wavenumber domain data Φ′(k _x , _ky , k _z , t=0) differs by a phase factor

Therefore, Φ'(k _x , _ky , k _z , t=0) can be obtained by interpolating Φ(k _x , _ky , z=0, f).

It should be emphasized here that the frequency domain data Φ(k _x , _ky , z=0, f) and the originally collected RF data are mutual Fourier transform pairs, while the spatial wavenumber domain data Φ′(k _x , k _y , k _z , t=0) and the image to be reconstructed are mutual Fourier transform pairs. A link between the known data and the image to be reconstructed is thus established. Preferably, the interpolation method is weighted inversion interpolation:

Φ'(k _x , _ky , k _z , t=0)=∫Φ(k _x , _ky , z=0, f)M(k _x , _ky , k _z , f)df, where:

Step (4), perform two-dimensional inverse Fourier transform on the spatial frequency domain data Φ'(k _x , _ky , k _z , t=0) to obtain spatial amplitude data s'(x, y, z, t =0).

Step (5), using the Hilbert demodulation method to demodulate the spatial amplitude data s' (x, y, z, t=0) to obtain envelope processing, and then perform logarithmic compression to adjust the contrast to obtain photoacoustic reconstruction image s*.

Step (6), calculate the tissue contrast CTR for s*, CTR is the ratio of the image brightness of the tissue in the image to the image brightness of the background position in the image, optimize the sound velocity distribution v(z'), and optimize the obtained V ₂ ( z') data are substituted into steps (3) to (5) to obtain a photoacoustic reconstructed image S ₂ * with the highest final tissue contrast CTR;

The specific steps for optimizing the sound velocity distribution v(z') are as follows:

(1), set the optimization times counter t=0, set the maximum optimization times T, the value range of T is 3-10, and randomly generate N sound velocity distributions V ₁ (z′) as the initial sample set K ₁ , the value of N is The value range is 300-1000, where the range of sound speed in each sound speed distribution is randomly within 1350-1650m/s;

(2) Reconstructing the photoacoustic image for each sound velocity distribution V ₁ (z′) in the above steps, and the reconstruction steps are the steps (1) to (5) in the above-mentioned photoacoustic tomography reconstruction method. According to the reconstruction method described above, N photoacoustic reconstructed images S ₁ * are obtained, and the CTR of each photoacoustic reconstructed image S ₁ * is calculated. The specific calculation method of CTR is the grayscale of the target object and the background in the photoacoustic reconstructed image S ₁ *. ratio of values;

(3) Sort the N calculated _CTRs in descending order of numerical value, select the sound velocity distribution corresponding to the first 30% larger CTRs to form the sample set Ka, and select the sound velocity distribution corresponding to the first 30%-60% larger CTRs to form The number of samples in the sample set K _b , K _a and K _b are all 0.3N;

(4) Randomly select a sound speed distribution sample V a (z') from Ka, randomly select a sound speed distribution sample V _{b (} z' _{) from K b ,} and take the mean value of the two sound speed distribution samples as the fused Sound velocity distribution V _c (z');

(5) Repeat step (4) to obtain 0.2N fused sound velocity distributions V _c (z′), and compare the unselected sound velocity distributions in step (4 _{) Ka and K b with} the fused sound velocity distribution V _c (z') together constitute a sample set K ₂ , and the number of samples in K ₂ is 0.4N;

(6), repeat steps (2)-(6) with K ₂ as the initial sample set, continuously update the sample set of the sound velocity distribution;

(7) When t=T, the optimization of the sound velocity distribution V(z′) is completed, and the sample with the largest CTR in the final sample set is selected as the final sound velocity distribution V ₂ (z′).

As shown in Figure 3, it is the reconstruction result of point sound source radio frequency data based on simulation. The original image is generated by k-wave simulation software, where the sampling frequency fs=30Mhz, the number of elements of the two-dimensional matrix acoustic transducer n= 192. The array element spacing of the two-dimensional matrix acoustic transducer is pt=0.2mm, and the speed of sound v=1540m/s. The radio frequency data R(x, y, z=0, t) is collected by the two-dimensional matrix acoustic transducer, and the time-reversal reconstruction and the reconstruction of the method of the present invention are carried out respectively; the results show that the lateral resolution of the reconstruction result of the method of the present invention is The rate is significantly better than the time-reversal reconstruction method.

Figure 4 shows a two-dimensional cross-section of the reconstruction result of the simulated blood vessel radio frequency data. The original data is generated by k-wave simulation software, where the sampling frequency is fs=30Mhz, and the number of elements of the two-dimensional matrix acoustic transducer is n. =192, the array element spacing of the two-dimensional matrix acoustic transducer is pt=0.2mm, and the speed of sound v=z+1500m/s. The radio frequency data R(x, y, z=0, t) is collected by the two-dimensional matrix acoustic transducer, and the filtered back-projection reconstruction and the reconstruction of the method of the present invention are respectively performed; the results show that the reconstructed image artifacts of the method of the present invention are far Much smaller than the filtered backprojection reconstruction, while the signal-to-noise ratio is significantly better than the filtered backprojection reconstruction method.

Figure 5 is a two-dimensional cross-section of the reconstruction result of the radio frequency data of the blood vessel model based on the real detection. The original data is actually collected by using the blood vessel model, where the sampling frequency is fs=50Mhz, and the number of elements of the two-dimensional matrix acoustic transducer is n=192 , The array element spacing of the two-dimensional matrix acoustic transducer is pt=0.25mm, and the sound speed is v=z+1500m/s. The radio frequency data R(x, y, z=0, t) is collected by the two-dimensional matrix acoustic transducer, and the filtered back-projection reconstruction and the reconstruction of the method of the present invention are respectively performed; the results show that the reconstructed image artifacts of the method of the present invention are far Much smaller than filtered back-projection reconstruction, while having better axial resolution.

It should be noted that the embodiments of the present invention and the features of the embodiments may be combined with each other under the condition of no conflict. Although this invention has been described to a certain extent, it will be apparent that suitable changes in various conditions may be made without departing from the spirit and scope of the invention. It is to be understood that the invention is not limited to the embodiments described, but is to be included within the scope of the claims, which include equivalents for each of the elements described.

Claims (2)

1. a photoacoustic tomography reconstruction method based on frequency domain wavenumber domain, is characterized in that, is made up of the following steps: (1) Acquire the radio frequency data R(x, y, z=0, t) collected by the two-dimensional matrix acoustic transducer, the system parameters of the two-dimensional matrix acoustic transducer and the sound velocity distribution V(z of the object to be measured) '), where the radio frequency data R(x, y, z=0, t) collected by the two-dimensional matrix acoustic transducer is collected in real time or saved in advance, and x and y are stored along the two-dimensional matrix acoustic transducer The position coordinates of x and y are perpendicular to each other, z is the position coordinates of the vertical two-dimensional matrix acoustic transducer direction, and t is the coordinates representing the time sequence; (2) Perform time gain compensation on the collected radio frequency data R(x, y, z=0, t) to obtain compensated radio frequency data S(x, y, z=0, t), where the time gain compensation method is: Multiply the value on each column of the radio frequency data R(x, y, z=0, t) by the weighting function Y(t), Y(t)=at+1, where the value range of a is 1-5; (3) Perform three-dimensional Fourier transform on the compensated radio frequency data S(x, y, z=0, t) along the x, y and t directions to obtain Φ(k _x , _ky , z=0, f), and then interpolate Φ(k _x , _ky , z=0, f) to obtain spatial wavenumber domain data Φ'(k _x , _ky , k _z , t=0), where the spatial wavenumber domain number Φ '(k _x , _ky , k _z , t=0)=∫Φ(kx, ky, z=0, f)M(k _x , _ky , k _z , f)df,

V(z') is the sound velocity distribution, k _x is the wave number along the x direction, _ky is the wave number along the y direction, k _z is the wave number along the z direction, and f is the frequency; (4), perform three-dimensional inverse Fourier transform on spatial frequency domain data Φ' (k _x , _ky , k _z , t=0) to obtain spatial amplitude data S' (x, y, z, t=0) ; (5) Demodulate the spatial amplitude data S' (x, y, z, t=0) to obtain envelope processing, and then perform contrast adjustment to obtain a photoacoustic reconstructed image S*, wherein the method of demodulating and obtaining the envelope It is the quadrature demodulation method or the Hilbert demodulation method, and the contrast adjustment method is the logarithmic compression method or the data obtained after de-envelope is made to the power of n, where 0<n<2; (6) Optimizing the sound velocity distribution V(z'), and substituting the optimized V ₂ (z') data into steps (3) to (5) to obtain a photoacoustic reconstructed image S with the highest final tissue contrast CTR ₂ *; 2. A kind of photoacoustic tomography reconstruction method based on frequency domain and wavenumber domain as claimed in claim 1, is characterized in that the concrete step of optimizing sound velocity distribution V(z') described in above-mentioned step (6) is: (1), set the optimization times counter t=0, set the maximum optimization times T, the value range of T is 3-10, and randomly generate N sound velocity distributions V ₁ (z′) as the initial sample set K ₁ , the value of N is The value range is 300-1000, where the range of sound speed in each sound speed distribution is randomly within 1350-1650m/s; (2) Reconstructing the photoacoustic image according to each sound velocity distribution V ₁ (z′) in the above steps, and the reconstruction step is the reconstruction method described in step (1) to step (5) in claim 1, to obtain N For the photoacoustic reconstructed image S ₁ *, the CTR of each photoacoustic reconstructed image S ₁ * is calculated, and the specific calculation method of the CTR is the ratio of the gray value of the target object and the background in the photoacoustic reconstructed image S ₁ *; (3) Sort the N calculated _CTRs in descending order of numerical value, select the sound velocity distribution corresponding to the first 30% larger CTRs to form the sample set Ka, and select the sound velocity distribution corresponding to the first 30%-60% larger CTRs to form The number of samples in the sample set K _b , K _a and K _b are all 0.3N; (4) Randomly select a sound speed distribution sample V a (z') from Ka, randomly select a sound speed distribution sample V _{b (} z' _{) from K b ,} and take the mean value of the two sound speed distribution samples as the fused Sound velocity distribution V _c (z'); (5) Repeat step (4) to obtain 0.2N fused sound velocity distributions V _c (z′), and compare the sound velocity distributions of unselected Ka and K _b in step (4 ₎ with the fused sound velocity distribution V _c (z') together constitute a sample set K ₂ , and the number of samples in K ₂ is 0.4N; (6), repeat steps (2)-(6) with K ₂ as the initial sample set, continuously update the sample set of the sound velocity distribution; (7) When t=T, the optimization of the sound velocity distribution V(z′) is completed, and the sample with the largest CTR in the final sample set is selected as the final sound velocity distribution V ₂ (z′).

Images