Object boundary extraction method suitable for ultrasonic and photoacoustic tomographyTechnical Field
The invention belongs to the technical field of ultrasonic tomography and photoacoustic tomography.
Background
The ultrasonic tomography is a novel ultrasonic imaging mode, which not only maintains the advantages of non-invasiveness, no radiation, good safety and the like of the traditional ultrasonic imaging, but also can provide echo images with high spatial resolution and quantitative images reflecting the sound velocity and attenuation of an imaging object due to the unique tomography imaging mode. Typical ultrasonic tomography systems use an annular ultrasonic transducer array to detect ultrasonic signals, emit scattered waves in the form of a synthetic aperture (SYNTHETIC APERTURE) to image an object, and the array elements of each channel sequentially emit ultrasonic waves and all channels together receive ultrasonic signals. The imaging mode has the advantages of more acquisition times, large data volume and low imaging speed. The plane wave emission mode can be used for scanning the whole imaging area in a short time, so that the acquisition time is greatly shortened, but the imaging quality is poor. Plane wave imaging can be achieved by using a polygonal linear array probe instead of an annular array, or by applying a delay to the annular array to cause it to emit plane waves. Photoacoustic tomography has the advantages of optical imaging and acoustic imaging, utilizes the light absorption difference of biological tissues, can provide structural and functional information of the inside of the tissues, and also frequently adopts an annular array to detect acoustic signals.
In the ultrasonic and photoacoustic tomography technology, boundary information of an imaging object is obtained a priori, so that the number of unknown variables in an image reconstruction area is reduced, the image reconstruction quality and speed are improved, and the method can be used for image post-processing and segmentation of the foreground and the background of an image. At present, boundary extraction in ultrasonic and photoacoustic tomography mainly utilizes digital image processing methods (such as Snakes, active contours, grabcut and other algorithms) to carry out boundary extraction on a reconstructed image, so that accurate boundary information cannot be obtained before image reconstruction, and further the image reconstruction process cannot be improved. Furthermore, the result of image boundary segmentation is severely dependent on the quality of the reconstructed image due to differences in imaging equipment and reconstruction parameters, and is susceptible to background noise and artifacts.
Disclosure of Invention
The invention aims to solve the technical problems that the boundary extraction method in ultrasonic and photoacoustic tomography cannot obtain accurate boundary information before image reconstruction, and the image boundary segmentation result is seriously dependent on the quality of reconstructed images and is easily influenced by background noise and artifacts.
In order to solve the technical problems, the technical scheme of the invention provides an object boundary extraction method suitable for ultrasonic and photoacoustic tomography, which is characterized by comprising the following steps:
Step 1, imaging an imaging target according to the requirements of ultrasonic tomography or photoacoustic tomography, receiving ultrasonic signals or photoacoustic signals through an ultrasonic transducer array, and selecting acquisition signals of a specific transducer receiving channel for subsequent boundary extraction according to imaging sequences generated by different imaging modes, wherein the acquisition signals comprise the following steps:
For ultrasonic tomography based on divergent wave emission, selecting proper emission-reception pair signals for subsequent boundary extraction;
For photoacoustic tomography, selecting acquisition signals of all transducer receiving channels for subsequent boundary extraction;
For ultrasonic tomography based on plane wave emission, plane waves of a certain angle are emitted, and acquisition signals of a transducer emission channel and a transducer receiving channel nearby the transducer emission channel are selected for subsequent boundary extraction;
Step 2, extracting the transit time of a first arrival echo from the acquired signals of the specific transducer receiving channel selected in the step 1;
step 3, let rrx＝(xrx,zrx be the spatial coordinate where a certain transducer receiving channel is located, r= (x, z) be the reflection point on the imaging target, c be the sound velocity of the coupling medium, and the transit time of the first arrival echo generated by the boundary is represented by the following formula (1), (2) or (3):
for ultrasound tomography based on divergent wave emission:
In the formula (1), rtx is the space coordinate of the divergent wave source, τSA(rtx,rrx, r) is the transit time of the first arrival echo generated by the boundary;
for photoacoustic tomography:
In equation (2), τPA(rrx, r) is the transit time of the first arrival echo generated by the boundary;
for ultrasound tomography based on plane wave emission:
in the formula (3), alpha is a plane wave emission angle, namely an included angle between a plane wave and the surface of a linear array ultrasonic probe, wherein the linear array ultrasonic probe is a real linear array ultrasonic probe or a virtual linear array ultrasonic probe obtained by overlapping and delaying other shaped ultrasonic arrays, L is the aperture size of the linear array ultrasonic probe, and tauPW(α,rrx, r) is the transit time of a first arrival echo generated by a boundary;
Wherein, the formulas (1), (2) or (3) are arranged into the form of an elliptic equation, representing all potential positions of the boundary points derived from the transit time of the first arrival echo, as shown in the following formula (4):
Ax2+Bxz+Cz2+Dx+Ez+F=0 (4)
In the formula (4), A-F is a constant;
Calculating a common tangent and corresponding tangent points for elliptic equations corresponding to two adjacent transducer receiving channels, and defining the tangent points as approximate boundary points and as initial boundary points;
Traversing the specific transducer receiving channels selected in the step 1 in sequence, and calculating to obtain all initial boundary points to obtain an initial boundary point set;
Step 4, post-processing is carried out on the initial boundary point set to obtain a final boundary detection result, and the method comprises the following steps:
step 401, distributing initial boundary points in an initial boundary point set into M fan-shaped areas which are uniformly distributed under a polar coordinate system according to a certain angle;
step 402, in each sector area, performing space outlier rejection on the initial boundary points, and only keeping the initial boundary point with the shortest distance from the polar coordinate origin;
Step 403, performing smooth denoising treatment on the discrete initial boundary points after outlier rejection;
and 404, interpolating the smoothed discrete initial boundary points, wherein all the obtained boundary points are final boundary detection results.
Preferably, in step 2, for any one segment of the acquired signal with the sampling length of N, the following steps are adopted to extract the transit time TOF of the first arrival echo of the acquired signal:
Step 201, calculating an AIC value of each sampling point, wherein the AIC value AIC (k) of the sampling point at the kth time is calculated according to the following formula (5):
in the formula (8), the amino acid sequence of the compound,AndThe variances of the two sections of acquired signals in the time periods [1, k ] and [ k+1, N ] are respectively;
In step 202, the transit time TOF of the first arrival echo is calculated based on the time corresponding to the sampling point with the minimum AIC value among all the sampling points, as shown in the following formula (6):
TOF=Ts*argmink{AIC(k)}+t0 (6)
In equation (9), Ts is a sampling interval, and T0 is a start time of an acquisition signal with a current sampling length N. The acquisition signals of the channels with lower signal-to-noise ratio can be discarded before detection, and the acquisition signals of the missing channels are interpolated by using the acquisition signals of other channels with high signal-to-noise ratio. In order to further reduce the interference of noise, the acquired signals may also be smoothed to improve the accuracy of subsequent steps.
Preferably, in step 3, the ultrasound tomography system based on divergent wave emission is composed of four modes:
The first mode) is that scattered waves are directly generated by a single array element or a plurality of array elements of the annular array ultrasonic transducer, wherein the annular array comprises a 360-degree closed annular array and also comprises an open annular array such as 90 degrees, and the description is omitted;
mode two), overlapping and delaying part of array elements of the annular ultrasonic array to obtain scattered waves generated by a virtual source;
The method comprises the steps of three) directly generating scattered waves by a single array element of the linear array ultrasonic transducer, or generating scattered waves by overlapping time delay of a plurality of array elements of the linear array ultrasonic transducer, and performing tomographic scanning by a rotary probe or a rotary imaging target;
mode IV), generating scattered waves by single or multiple array elements of the linear convex array ultrasonic transducer directly or through superposition time delay, and rotating a probe or rotating an imaging target to perform tomographic scanning.
Preferably, in step 3, the ultrasonic tomography system based on plane wave emission is composed of the following three modes:
the first mode is that a polygonal array is formed by a plurality of conventional linear array ultrasonic probes;
mode two), performing tomographic scanning by rotating or rotating an imaging target through a single or a plurality of conventional linear array ultrasonic probes;
Mode three) the annular ultrasonic array is overlapped and delayed, so that the annular ultrasonic array emits plane waves, and then:
And (3) calculating the formula under the coordinate system of the currently transmitted linear array ultrasonic probe or the virtual linear array ultrasonic probe obtained by overlapping and delaying the annular ultrasonic array.
Preferably, in step 3, for ultrasonic tomography based on plane wave emission, tangent point calculation is performed under the coordinate system where the emitted linear array probe is located, and then the coordinate is transformed to a uniform global coordinate system.
Preferably, in step 3, for ultrasound tomography based on divergent wave emission, A-F in equation (4) is calculated from rtx、rrx, c and τSA.
Preferably, in step 3, for ultrasonic tomography based on divergent wave emission, when the transducer receiving the signal is simultaneously the emission channel, then there is rtx＝rrx,A＝C＝1,B＝0,D＝-2xrx,E＝-2zrx,The formula (4) is simplified as:
In formula (7), τSA＝τSA(rtx＝rrx,rrx, r).
Preferably, in step 3, for photoacoustic tomography, a=c=1, b=0, d= -2xrx,E＝-2zrx,Equation (4) is expressed as:
in formula (8), τPA＝τPA(rrx, r).
Preferably, in step 3, for ultrasonic tomography based on plane wave emission, the A-F sum in formula (4) is calculated from α, L, rrx, c and τPW.
Preferably, in step 3, for ultrasound tomography based on plane wave emission, when plane wave angle a=0°, a=1, b=c=0, d= -2xrx,E＝2(cτPW-zrx,The formula (4) is simplified as follows;
In formula (9), τPW＝τPW(0°,rrx, r).
In step 4, preferably, when the space outlier is removed from the initial boundary points, the sector area is deflected according to the angle θ to obtain other sector areas, and then, the initial boundary points with the shortest distance from the origin point in all the sector areas after the angle deflection are reserved to remove the outlier.
Preferably, in step 4, in the post-processing process, parameters M and θ are adjusted according to the distribution of outliers, and a detection result closer to the boundary of the imaging target is retained.
The invention provides a method for directly extracting an imaging object boundary from an ultrasonic/photoacoustic channel signal, which is used for deducing the possible space position of the object boundary based on the transit time of a first arrival echo, accurately recovering the object boundary position by using methods such as space outlier rejection, smoothing, interpolation and the like, and is suitable for ultrasonic and photoacoustic tomography. Compared with a digital image processing method, the method provided by the invention is not influenced by the quality of the reconstructed image, and accurate boundary position information can be rapidly provided without reconstructing the image.
Drawings
FIG. 1 is a flow chart of a boundary extraction algorithm;
FIG. 2 is a schematic diagram of a coordinate system of ultrasonic tomography based on plane wave emission;
FIG. 3 is a schematic diagram of outlier rejection of boundary points;
FIG. 4 illustrates the transit time of a first arrival echo of finger ultrasound tomography;
FIG. 5 illustrates an ultrasonic tomography finger boundary detector based on scattered wave emission, wherein (a) an initial boundary point after tangent points are calculated for two adjacent receiving channels, (b) a boundary point after spatial outlier removal, (c) a discrete boundary point spatial position after outlier removal is smoothed, and (d) a final boundary detection result is obtained by interpolating the discrete boundary point spatial position;
fig. 6 illustrates the transit time of a first arrival echo of a finger photoacoustic tomography;
FIG. 7 illustrates photoacoustic tomography finger boundary detection, wherein (a) an initial boundary point after a tangent point is calculated for two adjacent receiving channels, (b) a final boundary detection result after post-processing;
Fig. 8 illustrates plane wave ultrasonic tomography breast boundary detection, in which (a) eight linear arrays form an octagonal ultrasonic tomography system, (b) a sound velocity distribution diagram of a breast model set in simulation, (c) an initial boundary point after tangent points are calculated for two adjacent receiving channels, and (d) a final boundary detection result after smooth interpolation processing.
Detailed Description
The application will be further illustrated with reference to specific examples. It is to be understood that these examples are illustrative of the present application and are not intended to limit the scope of the present application. Furthermore, it should be understood that various changes and modifications can be made by one skilled in the art after reading the teachings of the present application, and such equivalents are intended to fall within the scope of the application as defined in the appended claims.
As shown in fig. 1, the object boundary extraction method suitable for ultrasound and photoacoustic tomography provided by the invention comprises the following steps:
Step one, data acquisition
The imaging target is imaged according to the imaging requirements of ultrasound or photoacoustic, and ultrasound or photoacoustic signals are received through the ultrasound transducer array. According to different imaging sequences, selecting acquisition signals of a specific transducer receiving channel for subsequent boundary extraction, wherein:
For an ultrasonic tomography sequence based on divergent wave emission, selecting a proper emission and reception pair signal for subsequent boundary extraction, taking an acquisition signal of a receiving channel of a transducer which is both emission and reception as an example;
For the photoacoustic tomography sequence, selecting acquisition signals of all receiving channels for subsequent boundary extraction;
for an ultrasonic tomography sequence based on plane wave emission, plane waves (for example, 0° plane waves) with a certain angle are emitted, and acquisition signals of a transducer emission channel and a transducer receiving channel nearby the transducer emission channel are selected for subsequent boundary extraction.
Step two, extracting the transit time of the first arrival echo
The time-of-flight (TOF) of the first arrival echo is extracted from the acquired signals of the particular receive channel selected in step one using a red pool information content criterion (Akaike information criterion, AIC) or the like. For any segment of acquisition signal with a sampling length of N, there are:
the AIC value of each sampling point is calculated, wherein the AIC (k) value of the sampling point at the kth time is calculated as shown in the following formula (1):
In the formula (1), the components are as follows,AndThe variances of the two sections of acquired signals in the time periods [1, k ] and [ k+1, N ] are respectively;
and calculating the transit time TOF which is the first arrival echo based on the time corresponding to the sampling point with the minimum AIC value in all the sampling points, wherein the transit time TOF is shown in the following formula (2):
TOF=Ts*argmink{AIC(k)}+t0 (2)
In equation (2), Ts is a sampling interval, and T0 is a start time of an acquisition signal with a current sampling length N.
The transit time of the first arrival echo can also be calculated by adopting methods of first zero crossing point extraction, cross correlation, first echo peak value extraction and the like, and is mainly divided into the following four types:
1. Based on template matching, such as cross correlation, 2, based on evaluation functions, using different characteristics of the evaluation functions such as red pool information quantity criterion, energy ratio, entropy and the like near the transition time point as evaluation criteria, and 3, based on a neural network. 4. The method based on the signal characteristic points to be detected comprises the steps of first zero crossing point extraction, first echo peak value extraction, first significant peak value extraction and the like.
The detection result of the transit time is easily affected by the presence of noise and crosstalk in the original signal. The acquisition signals of the channels with lower signal-to-noise ratio can be discarded before detection, and the acquisition signals of the missing channels are interpolated by using the acquisition signals of other channels with high signal-to-noise ratio. In order to further reduce the interference of noise, the acquired signals may also be smoothed to improve the accuracy of subsequent steps.
Step three, initial boundary point extraction
A link of the boundary spatial position to the time of flight TOF of the first arriving echo is established. Assuming rrx＝(xrx,zrx) is the spatial coordinate where a certain receiving channel is located, r= (x, z) is the reflection point on the imaging object, c is the sound velocity of water, the transit time of the first arrival echo generated by the boundary can be expressed by the following equation (3), (4) or (5):
for ultrasound tomography based on divergent wave emission:
In the formula (3), rtx is the spatial coordinates of the receiving channel, τSA(rtx＝rrx,rrx, r) is the transit time of the first arrival echo generated by the boundary, wherein the ultrasound tomography system based on the emission of the scattered waves is composed of the following four modes:
Mode one) scattered waves are directly generated by single or multiple array elements of the annular array ultrasonic transducer, and note that;
mode two), overlapping and delaying part of array elements of the annular ultrasonic array to obtain scattered waves generated by a virtual source;
mode three), generating scattered waves directly by a single array element of the linear array ultrasonic transducer, or generating scattered waves by overlapping time delay of a plurality of array elements of the linear array ultrasonic transducer, and performing tomographic scanning by a rotary probe or a rotary imaging target;
mode IV), generating scattered waves by single or multiple array elements of the linear convex array ultrasonic transducer directly or through superposition time delay, and rotating a probe or rotating an imaging target to perform tomographic scanning. The expression (3) holds for all the above four modes.
For photoacoustic tomography:
In equation (4), τPA(rrx, r) is the transit time of the first arrival echo generated by the boundary;
for ultrasound tomography based on plane wave emission:
In equation (5), τPW(0°,rrx, r) is the transit time of the first arrival echo generated by the boundary, where a plane wave emission based ultrasound tomography system can be constructed in three ways:
the first mode is that a polygonal array is formed by a plurality of conventional linear array ultrasonic probes;
mode two), performing tomographic scanning by rotating or rotating an imaging target through a single or a plurality of conventional linear array ultrasonic probes;
mode three) the annular ultrasonic array is overlapped and delayed, so that the annular ultrasonic array emits plane waves.
Equation (5) holds for all three modes, and requires calculation under the coordinate system (local coordinate system) of the currently transmitted linear array ultrasonic probe (or the virtual linear array ultrasonic probe obtained by overlapping and delaying the annular ultrasonic array), as shown in fig. 2.
The finishing equations (3) - (5) are in the form of elliptic equations representing all potential positions of the boundary points derived from the transit time of the first arrival echo, wherein:
for ultrasound tomography based on divergent wave emission:
In formula (6), τSA＝τSA(rtx＝rrx,rrx, r);
for photoacoustic tomography:
In formula (7), τPA＝τPA(rrx, r);
for ultrasound tomography based on plane wave emission:
in formula (8), τPW＝τPW(0°,rrx, r).
And (3) calculating a common tangent and a corresponding tangent point for the elliptic equations corresponding to the two adjacent receiving channels as shown in the formulas (6), (7) or (8), and taking the tangent point as the approximation of the boundary point.
And (3) traversing the specific receiving channels selected in the step one in sequence, calculating tangential points, forming an initial boundary point set by all the calculated tangential points, and further defining each point in the initial boundary point set as an initial boundary point. For ultrasonic tomography based on plane wave emission, tangent point calculation needs to be carried out under the coordinate system (local coordinate system) where the emitted linear array probe is located, and after the tangent point calculation is completed, the coordinate is transformed to a unified global coordinate system.
Step four, post-treatment
And distributing the detected initial boundary points into M fan-shaped areas which are uniformly distributed at a certain angle under a polar coordinate system. And in each sector area, carrying out space outlier elimination on the initial boundary points, only keeping the initial boundary points with the shortest distance from the polar coordinate origin, wherein the number of the kept boundary points is determined by the number M of the set sector areas. In order to improve the outlier rejection effect, the angle of each sector area may be increased by θ without changing the number of sectors M. For example, in fig. 3, the black sector is angularly deflected by θ, resulting in two other sectors. And (3) reserving initial boundary points with shortest distances from the original points in all the fan-shaped areas after angle deflection, so as to realize outlier rejection.
And carrying out smooth denoising treatment on the discrete initial boundary points after outlier rejection.
And finally, interpolating the smoothed discrete initial boundary points, wherein all the obtained boundary points are final boundary detection results.
In the post-processing process, according to the distribution condition of outliers, the parameters M and theta can be adjusted to keep the detection result closer to the object boundary.
Example 1 ultrasonic tomography based on divergent wave emission
The object boundary extraction method suitable for ultrasonic and photoacoustic tomography disclosed in the embodiment comprises the following steps:
Step one, data acquisition
The finger was imaged ultrasonically using a 512 channel circular ultrasound array with a 40mm radius, a center frequency of 5MHz and a sampling frequency of 20MHz. Each channel array element sequentially excites scattered waves by using a synthetic aperture method, and all channels receive the scattered waves. An ultrasound echo image is reconstructed using 128 channels of data centered on the transmit array element, and finger boundaries are reconstructed using only the data received by the transmit array element.
Step two, extracting the transit time of the first arrival echo
The transit time of the first arrival echo is extracted using the AIC method, the result of which is shown in fig. 4.
Step three, initial boundary point extraction
According to formula (6), common tangents and tangent points are calculated for the elliptic equations between two adjacent receiving channels in sequence to obtain an initial set of boundary points, as shown in fig. 5 (a).
Step four, post-treatment
And performing space outlier rejection on the initial boundary point set, wherein the number M of the sector areas is set to 128, and the angle offset theta of the sector areas is set to 1 deg. As a result, as shown in fig. 5 (b), all outliers are eliminated, and the remaining boundary points are closer to the original boundary. In order to reduce the influence of noise, the discrete boundary points after outlier removal are subjected to spatial coordinate smoothing, and finally, a final boundary detection result is obtained through spatial coordinate interpolation, and the result is shown in fig. 5 (c) and (d).
Example 2 photoacoustic tomography
The object boundary extraction method suitable for ultrasonic and photoacoustic tomography disclosed in the embodiment comprises the following steps:
Step one, data acquisition
The finger was photoacoustic imaged using a 128 channel 30mm radius annular ultrasound array with a center frequency of 2.5MHz and a sampling rate of 40MHz. The laser emission pulse repetition frequency was 10Hz and the wavelength was 720nm. The 128-channel data are interpolated to 512 channels, so that the reconstruction quality of the photoacoustic image and the distribution density of initial boundary points are improved.
Step two, first arrival wave transit time extraction
The transit time of the first arriving echo is extracted using the AIC method and smoothed using a 7-point mean smoothing filter, the result of which is shown in fig. 6.
Step three, initial boundary point extraction
According to formula (7), common tangents and tangent points are calculated for each pair of elliptic equations between two adjacent receiving channels in turn, to obtain an initial set of boundary points, as shown in fig. 7 (a).
Step four, post-treatment
And (3) performing spatial outlier rejection on the initial boundary point set, and performing smoothing and interpolation processing on discrete coordinates after outlier rejection to obtain a final boundary detection result, as shown in fig. 7 (b).
Example 3 ultrasonic tomography based on plane wave emission
The object boundary extraction method suitable for ultrasonic and photoacoustic tomography disclosed in the embodiment comprises the following steps:
Step one, data acquisition
In simulation software k-Wave, a 1024-channel octagonal ultrasonic array with the radius of 40mm is used for carrying out plane Wave ultrasonic tomography on a numerical model of a breast slice, and the center frequency is set to be 4.5MHz. Eight linear arrays, each of which sequentially transmits 5 plane waves as shown in fig. 8 (a) and 8 (b), are wrapped around the breast model, the tilt angle of the transmitted plane waves being varied from-5 ° to 5 ° at intervals of 2.5 °, and three linear arrays centered on the transmit array receiving echo signals for image reconstruction as shown in fig. 2. And (4) performing boundary detection by using echo signals received by the transmitting array when the 0-degree plane wave is transmitted.
Step two, first arrival wave transit time extraction
The transit time of the first arriving echo is extracted using the AIC method and smoothed using a 5-point mean smoothing filter.
Step three, initial boundary point extraction
According to formula (8), common tangents and tangent points are calculated for each pair of elliptic equations between two adjacent receiving channels in turn, to obtain an initial set of boundary points, as shown in fig. 8 (c).
Step four, post-treatment
And (3) performing spatial outlier rejection on the initial boundary point set, and performing smoothing and interpolation processing on discrete coordinates after outlier rejection to obtain a final boundary detection result, as shown in fig. 8 (d).