CN117158911B - Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method - Google Patents
Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method Download PDFInfo
- Publication number
- CN117158911B CN117158911B CN202311387930.2A CN202311387930A CN117158911B CN 117158911 B CN117158911 B CN 117158911B CN 202311387930 A CN202311387930 A CN 202311387930A CN 117158911 B CN117158911 B CN 117158911B
- Authority
- CN
- China
- Prior art keywords
- sound
- ultrasonic
- ultrasonic transducer
- photoacoustic
- propagation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000003325 tomography Methods 0.000 title claims abstract description 47
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 22
- 230000031700 light absorption Effects 0.000 claims abstract description 16
- 238000003384 imaging method Methods 0.000 claims abstract description 10
- 230000001902 propagating effect Effects 0.000 claims abstract description 7
- 239000013598 vector Substances 0.000 claims description 30
- 238000010521 absorption reaction Methods 0.000 claims description 28
- 238000002604 ultrasonography Methods 0.000 claims description 22
- 238000004364 calculation method Methods 0.000 claims description 16
- 239000000523 sample Substances 0.000 claims description 16
- 238000010895 photoacoustic effect Methods 0.000 claims description 10
- 230000003044 adaptive effect Effects 0.000 claims description 8
- 238000005452 bending Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000000644 propagated effect Effects 0.000 claims description 3
- 210000001519 tissue Anatomy 0.000 description 51
- XLYOFNOQVPJJNP-ZSJDYOACSA-N Heavy water Chemical compound [2H]O[2H] XLYOFNOQVPJJNP-ZSJDYOACSA-N 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 210000004209 hair Anatomy 0.000 description 4
- 230000006978 adaptation Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- INGWEZCOABYORO-UHFFFAOYSA-N 2-(furan-2-yl)-7-methyl-1h-1,8-naphthyridin-4-one Chemical compound N=1C2=NC(C)=CC=C2C(O)=CC=1C1=CC=CO1 INGWEZCOABYORO-UHFFFAOYSA-N 0.000 description 1
- 206010006187 Breast cancer Diseases 0.000 description 1
- 208000026310 Breast neoplasm Diseases 0.000 description 1
- 108010064719 Oxyhemoglobins Proteins 0.000 description 1
- 239000006096 absorbing agent Substances 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 210000004204 blood vessel Anatomy 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 108010002255 deoxyhemoglobin Proteins 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000001727 in vivo Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007102 metabolic function Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012634 optical imaging Methods 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Landscapes
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
The invention belongs to the technical field of photoacoustic image reconstruction, and discloses a multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method, which comprises the following steps that firstly, ultrasonic transducers are used for transmitting ultrasonic waves to irradiate biological tissues in an ultrasonic tomography equipment system; step two, the ultrasonic transducer receives ultrasonic signals sent by biological tissues; thirdly, reconstructing sound velocity distribution of the biological tissue by adopting an iterative reconstruction algorithm of a curved ray model; dividing biological tissue into a plurality of areas with different sound speeds; step five, emitting pulse laser light to irradiate the tissue by a laser in the photoacoustic tomography equipment system; step six, the ultrasonic transducer receives the photoacoustic signal sent by the light absorption point in the tissue; step seven, calculating the multi-sound speed involved in the process of propagating the photoacoustic signal to the ultrasonic transducer; and step eight, the multi-sound velocity self-adaption is applied to a delay summation algorithm, and a photoacoustic tomography image is reconstructed. The invention effectively improves the quality of the photoacoustic imaging reconstructed image.
Description
Technical Field
The invention belongs to the technical field of photoacoustic image reconstruction, and particularly relates to a multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction system and method.
Background
The photoacoustic tomography is a nondestructive medical imaging method developed in recent years, combines the high-contrast characteristic of pure optical imaging and the high-penetration depth characteristic of pure ultrasonic imaging, can provide living tissue imaging with high resolution and high contrast, provides an important means for researching the structural morphology, physiological characteristics, metabolic functions, case characteristics and the like of biological tissues, and has wide application prospects in the fields of biomedical clinical diagnosis and in-vivo tissue structure and functional imaging. The main principle of photoacoustic tomography is the conversion of light energy into acoustic energy, the photoacoustic effect. When a beam of pulsed laser irradiates biological tissue, the light absorber inside the tissue absorbs the laser energy, and instantaneous thermal elastic expansion occurs to excite ultrasonic waves, namely photoacoustic signals. The photoacoustic signal propagates toward the tissue surface and is ultimately received by the ultrasound transducer. According to the time delay and amplitude of the received photoacoustic signal and the preset sound velocity of the system, the optical characteristic distribution, particularly the light absorption coefficient distribution, inside biological tissues is reconstructed after the computer processing. These information items can be used for quantitative measurement of specific substances in biological tissues, such as oxyhemoglobin and deoxyhemoglobin in blood, and on the basis of this, clinical studies such as breast cancer and the like are carried out.
In view of complexity of biological tissues, in order to simplify the problem in photoacoustic tomography, existing methods generally assume that propagation velocity of a photoacoustic signal in biological tissues is constant, ignoring influence caused by uneven sound velocity distribution in tissues. Studies have shown that biological tissues with different compositions often have different acoustic properties, where the difference in sound velocity inside soft tissues can reach 10%. Therefore, inaccuracy of the preset sound velocity in the photoacoustic tomography image reconstruction algorithm is one of the main causes of poor quality of the reconstructed image.
Disclosure of Invention
In order to solve the technical problems, the specific technical scheme of the multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method is as follows:
a multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method comprises the following steps:
step one, an ultrasonic transducer emits ultrasonic waves to irradiate biological tissues;
step two, the ultrasonic transducer receives ultrasonic signals sent by biological tissues;
thirdly, reconstructing sound velocity distribution of the biological tissue by adopting an iterative reconstruction algorithm of a curved ray model;
dividing biological tissue into a plurality of areas with different sound speeds;
step five, the laser emits pulse laser to irradiate the tissue;
step six, the ultrasonic transducer receives the photoacoustic signal sent by the light absorption point in the tissue;
step seven, calculating the multi-sound speed involved in the process of propagating the photoacoustic signal to the ultrasonic transducer;
and step eight, the multi-sound velocity self-adaption is applied to a delay summation algorithm, and a photoacoustic tomography image is reconstructed.
Further, in the first step, the ultrasonic transducer receives ultrasonic waves and emits ultrasonic waves, the ultrasonic transducer array is annular, and a plurality of ultrasonic transducers are uniformly selected from the ultrasonic transducer array to emit ultrasonic waves successively.
Further, the ultrasonic transducer in the second step sequentially receives the ultrasonic waves emitted by the biological tissue to be tested in the first step, the ultrasonic waves reach the ultrasonic transducer through the acoustic wave lens in the ultrasonic probe, and the ultrasonic transducer records the sound pressure.
Further, the iterative reconstruction algorithm of the curved ray model in the third step is as follows:
the bending rays are used for representing an array element in the annular ultrasonic transducer array to emit ultrasonic signals, the ultrasonic waves are refracted when passing through different sound speed areas, and finally the ultrasonic transducer receives the refraction process, and the formula is as follows:
，
wherein,a first arrival time vector for the ultrasonic signal to reach the ultrasonic transducer; />Is the slowness vector of the biological tissue to be tested; />To link the coefficient matrix of the two with the spatial distribution of the ultrasound transducer and the slowness distribution of the tissue +.>In the process of sound velocity reconstruction, the theoretical sound wave arrival time and the experimental observation arrival time based on the bending ray model are continuously compared until the deviation of the theoretical sound wave arrival time and the experimental observation arrival time is smaller than a set threshold value.
Further, the fourth step includes the following specific steps:
and (3) traversing each sound source point according to the sound velocity distribution reconstructed in the step (III), classifying other sound source points which are different from the sound velocity of the point within a set range into the same region, traversing each region, taking the average sound velocity of all sound source points in the region as the sound velocity of the region, and dividing a biological tissue into a plurality of regions with different sound velocities.
Further, the irradiated biological tissue in the fifth step generates a photoacoustic effect, which satisfies the equation:
， (1)
wherein the method comprises the steps ofAcoustic pressure generated for thermal expansion of biological tissue, < >>Is specific heat capacity->Is an isobaric volume thermal expansion coefficient, +.>Is the empirical sound velocity,/->For the heat source function excited by the incident pulsed laser, +.>，/>Is a spatial light absorption function of biological tissue, +.>For the light irradiation function, the laser strikes the absorption point and the ultrasonic transducer starts to receive signals at the same time.
Further, the pulse laser irradiates the absorption point below or beside the biological tissue to be measured, and the plane of the ultrasonic transducer array is perpendicular to the plane of the laser.
Further, the ultrasonic transducer adopted in the step six is consistent with the ultrasonic transducer adopted in the step two, the ultrasonic transducer receives the ultrasonic waves generated by the photoacoustic effect in the step five, the ultrasonic waves reach the ultrasonic transducer through the acoustic lens in the ultrasonic probe, the sound pressure intensity is recorded, and the number of the ultrasonic transducers is set asRecord->The sound pressure signal of each ultrasonic transducer generated by the photoacoustic effect in equation (1) is +.>, />。
Further, the photoacoustic signal in the seventh step passes through a plurality of regions with different sound speeds in the process of propagating from the light absorption point to the tissue surface, and is finally received by the ultrasonic transducer, on the basis of a plurality of regions with different sound speeds marked in the fourth step, the sound speed region of the tissue is discretized into grids under the limitation of the imaging region, the photoacoustic signal is propagated to the ultrasonic transducer along a straight line, and the propagation vector is marked asWherein->Indicating the light absorption point +.>Representing an ultrasound transducer, while->Is the origin of the propagation path, < > and>is the end point of the propagation path, and when the propagation path of the photoacoustic signal is calculated in each region, the propagation vector is selected +.>And calculating each propagation path and the sound velocity corresponding to the propagation path by the path closest to the same, and finally obtaining multi-sound velocity information on the whole propagation path.
Further, the step eight applies multi-sound speed self-adaption to the calculation of the propagation time of the photoacoustic signal from the absorption point to the ultrasonic transducer, and for each absorption point, a delay summation algorithm formula is as follows:
，
wherein,representing a photoacoustic tomography image reconstructed based on a delay-and-sum algorithm, < >>Representing ultrasound transducer->At time->Received photoacoustic signal,/->Representing the propagation of a photoacoustic signal from the absorption point to the ultrasound transducer>The desired propagation time is set to be the desired propagation time,is the weight of the sample, and the weight of the sample,
adaptive application of multiple sound speeds to propagation timeIn the calculation of (2), the calculation formula is as follows:
，
wherein the method comprises the steps ofIs the required propagation time,/->For the absorption point, < > for>Is an ultrasonic transducer->At sound speed +.>The propagation path over the region of (2) is->，/>Is the boundary point of the region, vector +.>Nearest propagation vector>I.e.And->The included angle between the two vectors is the smallest; similarly, at sound speed +.>The propagation path over the region of (2) is->，/>Is the boundary point of the region, likewise vector +.>Nearest propagation vector>The method comprises the steps of carrying out a first treatment on the surface of the Up to the end of the propagation path->A region where the sound velocity is +.>Propagation path is->The method comprises the steps of carrying out a first treatment on the surface of the The weight formula is:
，
wherein,is the distance of the absorption point to the ultrasound transducer, +.>Is the included angle between the absorption point and the ultrasonic transducer, in the two-dimensional X-Y plane, the circle center of the annular ultrasonic transducer array is taken as the origin of coordinates, and the absorption point is assumed to be +.>Coordinates areUltrasound transducer->Coordinates of->Then->The calculation formula is as follows:
，
the calculation formula is as follows:
。
the multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method has the following advantages:
according to the invention, the sound velocity distribution of biological tissues is reconstructed by utilizing ultrasonic sound velocity tomography, different sound velocity areas are divided, the multi-sound velocity involved in the propagation process of the photoacoustic signal is calculated on the basis, and the multi-sound velocity is adaptively applied to the reconstruction of the photoacoustic tomography image, so that the adverse effect of inaccuracy of the preset sound velocity on the photoacoustic imaging is avoided, and the quality of the photoacoustic imaging reconstructed image is effectively improved.
Drawings
FIG. 1 is a flow chart of the method of the present invention;
FIG. 2 is a flow chart of an iterative reconstruction algorithm based on a curved ray model in the method of the present invention;
FIG. 3 is a schematic diagram of calculating propagation paths of photoacoustic signals in different acoustic velocity regions in the method of the present invention;
FIG. 4 is a graph of experimental samples for simulation verification using k-wave;
fig. 5 (a) is a view of the effect of photoacoustic tomography reconstruction using a preset sound speed;
fig. 5 (b) is a view showing the effect of reconstruction of a multi-sound-velocity adaptive photoacoustic tomography image using the present invention.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
Considering that the photoacoustic tomography system and the ultrasonic tomography system have similar acoustic signal receiving and processing procedures, the fusion of the two systems can be realized with lower software and hardware cost. Therefore, in order to solve the problem of inaccurate preset sound velocity in the existing photoacoustic tomography technology, the invention provides a photoacoustic tomography image reconstruction algorithm which does not depend on the preset sound velocity and can be self-adaptive to different biological tissues based on ultrasonic sound velocity tomography, the algorithm firstly reconstructs sound velocity distribution of the biological tissues based on ultrasonic sound velocity tomography, different sound velocity areas are divided, the multi-sound velocity involved in the propagation path from a photoacoustic signal to an ultrasonic transducer (an ultrasonic sensor) is calculated on the basis, and finally the multi-sound velocity is self-adaptive to the photoacoustic tomography image reconstruction, so that an image with higher quality is reconstructed.
As shown in fig. 1, the invention discloses a multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method, which comprises the following steps:
step one, transmitting ultrasonic waves to irradiate biological tissues by using an ultrasonic transducer in an ultrasonic tomography equipment system;
the ultrasonic tomography equipment system can complete the reconstruction of the photoacoustic tomography image and the reconstruction of the ultrasonic tomography image, and an ultrasonic transducer used by the system can receive ultrasonic waves and emit ultrasonic waves. The ultrasonic transducer array in the system is annular, surrounds the biological tissue to be tested, and is filled with heavy water. And uniformly selecting a plurality of ultrasonic transducers from the ultrasonic transducer array to sequentially emit ultrasonic waves.
Step two, the ultrasonic transducer receives ultrasonic signals sent by biological tissues;
the ultrasonic transducer successively receives ultrasonic waves emitted by the biological tissue to be tested in the first step, the ultrasonic waves reach the ultrasonic transducer through the acoustic wave lens in the ultrasonic probe, and the sound pressure is recorded. The greater the number of ultrasound transducers used, the higher the quality of the reconstructed image.
Thirdly, reconstructing sound velocity distribution of the biological tissue by adopting an iterative reconstruction algorithm of a curved ray model;
the algorithm principle is shown in fig. 2, one array element in the annular ultrasonic transducer array emits an ultrasonic signal, the ultrasonic wave can be refracted when passing through different sound speed areas, and finally the ultrasonic wave is received by the ultrasonic transducer, and the refraction process can be characterized by using bending rays, wherein the following formula is shown as follows:
，
wherein,a first arrival time vector for the ultrasonic signal to reach the ultrasonic transducer; />Is the slowness (reciprocal of sound velocity) vector of the measured biological tissue; />To link the coefficient matrix of the two with the spatial distribution of the ultrasound transducer and the slowness distribution of the tissue +.>In the process of sound velocity reconstruction, the theoretical sound wave arrival time and the experimental observation arrival time based on the bending ray model are continuously compared until the deviation of the theoretical sound wave arrival time and the experimental observation arrival time is smaller than a set threshold value.
Dividing biological tissue into a plurality of areas with different sound speeds;
and (3) traversing each sound source point according to the sound velocity distribution reconstructed in the step three, and classifying other sound source points which are different from the sound velocity of the point within a certain range (such as < = 10 m/s) into the same region. Finally, each region is traversed, and the average sound velocity of all sound source points in the region is taken as the sound velocity of the region. Further, one biological tissue is divided into a plurality of regions having different sound speeds.
Step five, emitting pulse laser to irradiate the tissue in an ultrasonic tomography equipment system;
in the equipment system of ultrasonic tomography, the irradiated biological tissue generates a photoacoustic effect, and the equation is satisfied:
， (1)
wherein the method comprises the steps ofAcoustic pressure generated for thermal expansion of biological tissue, < >>Is specific heat capacity->Is an isobaric volume thermal expansion coefficient, +.>Is the empirical sound velocity,/->For the heat source function excited by the incident pulsed laser, +.>，/>Is a spatial light absorption function of biological tissue, +.>For the light irradiation function, the laser strikes the absorption point and the ultrasonic transducer starts to receive signals at the same time.
In the equipment system of ultrasonic tomography, the working wavelength of the laser pulse emitted by the laser is 1064nm,maximum power of about 2J/cm 2 . The pulsed laser irradiates the absorption point below or beside the biological tissue to be measured. The ultrasonic transducer surrounds the biological tissue under test, filling heavy water between them. The plane in which the ultrasound transducer array is located is generally perpendicular to the plane in which the laser is located.
Step six, the ultrasonic transducer receives the photoacoustic signal sent by the light absorption point in the tissue;
and step six, the ultrasonic transducer adopted in the step two is consistent with the ultrasonic transducer adopted in the step two, and heavy water is filled between the ultrasonic transducer and the biological tissue around the detected biological tissue. The ultrasonic transducer receives the ultrasonic wave generated by the photoacoustic effect in the fifth step, the ultrasonic wave reaches the ultrasonic transducer through the acoustic wave lens in the ultrasonic probe, and the sound pressure intensity is recorded. The greater the number of ultrasound transducers used, the higher the quality of the reconstructed image. The ultrasonic transducers are arranged in a ring shape. Let the number of ultrasonic transducers beRecord->The sound pressure signal of each ultrasonic transducer generated by the photoacoustic effect in equation (1) is +.>, />。
Step seven, calculating the multi-sound speed involved in the process of propagating the photoacoustic signal to the ultrasonic transducer;
the photoacoustic signal in the seventh step passes through a plurality of areas with different sound speeds in the process of propagating from the light absorption point to the tissue surface and is finally received by the ultrasonic transducer. On the basis of a plurality of different sound speed areas marked in the fourth step, the sound speed area of the tissue is limited by the imaging area and is discretized into a grid with a certain size, the photoacoustic signal is assumed to be propagated to the ultrasonic transducer along a straight line, and the propagation vector is recorded asWherein->Indicating the light absorption point +.>Representing an ultrasound transducer, while->Is the origin of the propagation path, < > and>is the end point of the propagation path, and when the propagation path of the photoacoustic signal is calculated in each region, the propagation vector is selected +.>The path closest to the coincidence. As each propagation path is calculated, the sound velocity corresponding to the propagation path is also calculated, and finally multi-sound velocity information on the entire propagation path is obtained.
The propagation path of the photoacoustic signal in each sound velocity region needs to be calculated separately, as shown in fig. 3, from the light absorption pointPropagate to ultrasonic transducer->Propagation vector is +.>The propagation path of the photoacoustic signal in the sound velocity region 2 is calculated as +>This is because of the vector +.>And propagation vector->The included angle of (c) is the smallest, i.e. closest to the propagation direction. Similarly, calculateThe propagation path in the sound velocity region 3 is +.>The propagation path in region 1 is +.>. Finally, the total propagation path of the photoacoustic signal is +.>The corresponding multi-sound speed is +.>。
Step eight, the multi-sound speed self-adaption is applied to a delay summation algorithm, and a photoacoustic tomography image is reconstructed;
the multi-sound speed adaptation is applied to the calculation of the propagation time of a photoacoustic signal from an absorption point (sound source point) to an ultrasonic transducer. For each absorption point, the delay-and-sum algorithm formula is as follows:
，
wherein,representing a photoacoustic tomography image reconstructed based on a delay-and-sum algorithm, < >>Representing ultrasound transducer->At time->Received photoacoustic signal,/->Representing the propagation of a photoacoustic signal from the absorption point to the ultrasound transducer>The desired propagation time is set to be the desired propagation time,is a weight.
Adaptive application of multiple sound speeds to propagation timeIn the calculation of (2), the calculation formula is as follows:
，
wherein the method comprises the steps ofIs the required propagation time,/->For the absorption point, < > for>Is an ultrasonic transducer->At sound speed +.>The propagation path over the region of (2) is->，/>Is the boundary point of the region, vector +.>Nearest propagation vector>I.e. +.>And->The included angle between the two vectors is the smallest; similarly, at sound speed +.>The propagation path over the region of (2) is->，/>Is the boundary point of the region, likewise vector +.>Nearest propagation vector>The method comprises the steps of carrying out a first treatment on the surface of the Up to the end of the propagation path->A region where the sound velocity is +.>Propagation path is->. Considering that the signal intensity received by the ultrasonic transducer and the included angle between the absorption point and the ultrasonic transducer form a cosine relation, and meanwhile, the signal intensity is inversely proportional to the distance, the weight formula in the invention is designed as follows:
，
wherein,is the distance of the absorption point to the ultrasound transducer, +.>Is the included angle between the absorption point and the ultrasonic transducer, and in the two-dimensional X-Y plane, the annular ultrasonic transducer is usedThe center of the array is the origin of coordinates, assuming the absorption point +.>Coordinates areUltrasound transducer->Coordinates of->Then->The calculation formula is as follows:
，
as can be seen from the calculation formula of the included angle between the two vectors in the analysis geometry,the calculation formula is as follows:
。
the invention will be further described with reference to a simulated application example.
The simulation examples are as follows:
the experimental sample shown in fig. 4 is generated by adopting a k-wave simulation tool, three hairs (marked by 1, 2 and 3 in the figure) buried in the sample and formed into a blood vessel pattern form an acoustic source point and also are light absorption points, the sample is divided into three areas with different sound speeds, the sound speed of the area 1 is 1400m/s, the sound speed of the area 2 is 1510m/s, the sound speed of the area 3 is 1590m/s, and the three hairs are completely contained in the area 3. And generating a semi-annular ultrasonic transducer array by adopting a k-wave simulation tool, wherein the semi-annular ultrasonic transducer array comprises 512 ultrasonic transducers which are uniformly distributed on the semi-ring.
And similarly, simulating ultrasonic emission by adopting a k-wave simulation tool, receiving ultrasonic signals by an ultrasonic transducer array, and sending the ultrasonic signals into an upper computer to reconstruct sound velocity distribution of a sample and divide a sound velocity region.
Similarly, after a k-wave simulation tool is used for simulating a pulse laser to irradiate a sample, the hair emits a photoacoustic signal and propagates outwards, an ultrasonic transducer array receives the photoacoustic signal and sends the photoacoustic signal into an upper computer for reconstructing a multi-sound-velocity self-adaptive photoacoustic tomography image, and the reconstructed image is shown in fig. 5 (a).
For comparison of photoacoustic imaging effects of the method proposed by the present invention and the method of presetting sound velocity, the preset sound velocity is set to 1500m/s (the average sound velocity of three sound velocity regions in the above sample), and the same is substituted into a delay-and-sum algorithm for reconstruction of a photoacoustic tomography image, the reconstructed image being as shown in fig. 5 (b).
As can be seen from fig. 5 (a) and fig. 5 (b), compared with the method of presetting the sound velocity, the photoacoustic tomography image obtained by the method has obviously reduced artifacts and distortion, three hairs in the image are clearly visible, and the image quality is obviously improved.
The above description of the embodiments is only for aiding in the understanding of the method of the present invention and its core concept, and it should be noted that it will be apparent to those skilled in the art that several modifications and adaptations can be made without departing from the scope of the present invention.
Claims (6)
1. The multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method is characterized by comprising the following steps of:
step one, an ultrasonic transducer emits ultrasonic waves to irradiate biological tissues;
the ultrasonic transducer is used for receiving ultrasonic waves and transmitting the ultrasonic waves, the ultrasonic transducer array is annular, and a plurality of ultrasonic transducers are uniformly selected from the ultrasonic transducer array to sequentially transmit the ultrasonic waves;
step two, the ultrasonic transducer receives ultrasonic signals sent by biological tissues;
thirdly, reconstructing sound velocity distribution of the biological tissue by adopting an iterative reconstruction algorithm of a curved ray model;
dividing biological tissue into a plurality of areas with different sound speeds;
aiming at the sound velocity distribution reconstructed in the step three, traversing each sound source point, classifying other sound source points which have the sound velocity difference within a set range from the sound velocity of the point into the same region, traversing each region, taking the average sound velocity of all sound source points in the region as the sound velocity of the region, and dividing a biological tissue into a plurality of regions with different sound velocities;
step five, the laser emits pulse laser to irradiate the tissue;
step six, the ultrasonic transducer receives the photoacoustic signal sent by the light absorption point in the tissue;
step seven, calculating the multi-sound speed involved in the process of propagating the photoacoustic signal to the ultrasonic transducer;
the photoacoustic signal passes through a plurality of areas with different sound speeds in the process of propagating from the light absorption point to the tissue surface and is finally received by the ultrasonic transducer, on the basis of a plurality of areas with different sound speeds marked in the fourth step, the sound speed area of the tissue is scattered into grids under the limitation of the imaging area, the photoacoustic signal is propagated to the ultrasonic transducer along a straight line, and the propagation vector is recorded asWherein->Indicating the light absorption point +.>Representing an ultrasound transducer, while->Is the origin of the propagation path, < > and>is the end point of the propagation path, and when the propagation path of the photoacoustic signal is calculated in each region, the propagation vector is selected +.>Calculating each propagation path and sound velocity corresponding to the propagation path by the path closest to the same, and finally obtaining multi-sound velocity information on the whole propagation path;
step eight, the multi-sound speed self-adaption is applied to a delay summation algorithm, and a photoacoustic tomography image is reconstructed;
applying multi-sound speed self-adaption to propagation time calculation of a photoacoustic signal from an absorption point to an ultrasonic transducer, a delay summation algorithm formula for each absorption point is as follows:
，
wherein,representing a photoacoustic tomography image reconstructed based on a delay-and-sum algorithm, < >>Representing ultrasound transducer->At time->Received photoacoustic signal,/->Representing the propagation of a photoacoustic signal from the absorption point to the ultrasound transducer>The desired propagation time is set to be the desired propagation time,is the weight of the sample, and the weight of the sample,
adaptive application of multiple sound speeds to propagation timeIn the calculation of (2), the calculation formula is as follows:
，
wherein the method comprises the steps ofIs the required propagation time,/->For the absorption point, < > for>Is an ultrasonic transducer->At sound speed +.>The propagation path over the region of (2) is->，/>Is the boundary point of the region, vector +.>Nearest propagation vector>I.e. +.>Andthe included angle between the two vectors is the smallest; similarly, at the speed of soundIs->The propagation path over the region of (2) is->，/>Is the boundary point of the region, likewise vector +.>Nearest propagation vector>The method comprises the steps of carrying out a first treatment on the surface of the Up to the end of the propagation path->A region where the sound velocity is +.>Propagation path is->The method comprises the steps of carrying out a first treatment on the surface of the The weight formula is:
，
wherein,is the distance of the absorption point to the ultrasound transducer, +.>Is the included angle between the absorption point and the ultrasonic transducer, in the two-dimensional X-Y plane, the circle center of the annular ultrasonic transducer array is taken as the origin of coordinates, and the absorption point is assumed to be +.>Coordinates areUltrasound transducer->Coordinates of->Then->The calculation formula is as follows:
，
the calculation formula is as follows:
。
2. the method for reconstructing a multi-sound-velocity adaptive photoacoustic tomography image according to claim 1, wherein the ultrasonic transducer in the second step sequentially receives the ultrasonic waves emitted from the biological tissue to be measured in the first step, the ultrasonic waves reach the ultrasonic transducer through the acoustic lens in the ultrasonic probe, and the ultrasonic transducer records the sound pressure.
3. The multi-sound-velocity adaptive photoacoustic tomography image reconstruction method according to claim 1, wherein the iterative reconstruction algorithm of the curved ray model in the third step is as follows:
the bending rays are used for representing an array element in the annular ultrasonic transducer array to emit ultrasonic signals, the ultrasonic waves are refracted when passing through different sound speed areas, and finally the ultrasonic transducer receives the refraction process, and the formula is as follows:
，
wherein,a first arrival time vector for the ultrasonic signal to reach the ultrasonic transducer; />Is the slowness vector of the biological tissue to be tested; />To link the coefficient matrix of the two with the spatial distribution of the ultrasound transducer and the slowness distribution of the tissue +.>In the process of sound velocity reconstruction, the theoretical sound wave arrival time and the experimental observation arrival time based on the bending ray model are continuously compared until the deviation of the theoretical sound wave arrival time and the experimental observation arrival time is smaller than a set threshold value.
4. The method of claim 1, wherein the irradiated biological tissue in the fifth step generates a photoacoustic effect, satisfying the equation:
， (1)
wherein the method comprises the steps ofAcoustic pressure generated for thermal expansion of biological tissue, < >>Is specific heat capacity->Is an isobaric volume thermal expansion coefficient, +.>Is the empirical sound velocity,/->For the heat source function excited by the incident pulsed laser, +.>，/>Is a spatial light absorption function of biological tissue, +.>For the light irradiation function, the laser strikes the absorption point and the ultrasonic transducer starts to receive signals at the same time.
5. The method for reconstructing a multisonic adaptive photoacoustic tomography image of claim 4, wherein the pulsed laser irradiates the absorption points below or beside the biological tissue to be measured, and the plane in which the ultrasonic transducer array is located is perpendicular to the plane in which the laser is located.
6. The method for reconstructing a multi-sound-velocity adaptive photoacoustic tomography image according to claim 4, wherein the ultrasonic transducer used in the sixth step is identical to the ultrasonic transducer used in the second step, the ultrasonic transducer receives the ultrasonic waves generated by the photoacoustic effect in the fifth step, and the ultrasonic waves reach the ultrasonic transducer through the acoustic lens in the ultrasonic probe, and the sound pressure intensity is recorded, and the number of the ultrasonic transducers is set to beRecord->Individual ultrasonic transducersThe sound pressure signal generated by the photoacoustic effect in equation (1) is +.>, 。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311387930.2A CN117158911B (en) | 2023-10-25 | 2023-10-25 | Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311387930.2A CN117158911B (en) | 2023-10-25 | 2023-10-25 | Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117158911A CN117158911A (en) | 2023-12-05 |
| CN117158911B true CN117158911B (en) | 2024-01-23 |
Family
ID=88943395
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311387930.2A Active CN117158911B (en) | 2023-10-25 | 2023-10-25 | Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117158911B (en) |
Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101214156A (en) * | 2008-01-10 | 2008-07-09 | 复旦大学 | Acoustic velocity inhomogeneous medium thermoacoustic imaging reconstruction algorithm |
| JP2012075464A (en) * | 2010-09-30 | 2012-04-19 | Fujifilm Corp | Photoacoustic diagnostic imaging equipment, image generating method, and program |
| CN103445765A (en) * | 2013-09-24 | 2013-12-18 | 南京大学 | Acoustic velocity correction method for photoacoustic imaging |
| CN104545811A (en) * | 2014-12-26 | 2015-04-29 | 深圳先进技术研究院 | Intravascular imaging system and method |
| WO2015175431A1 (en) * | 2014-05-12 | 2015-11-19 | University Of Washington | Real-time photoacoustic and ultrasound imaging system and method |
| WO2015189268A2 (en) * | 2014-06-10 | 2015-12-17 | Ithera Medical Gmbh | Device and method for hybrid optoacoustic tomography and ultrasonography |
| CN105249993A (en) * | 2015-11-16 | 2016-01-20 | 南京大学 | Method for selecting optimum sound velocity group to optimize ultrasonic imaging through photoacoustic imaging |
| JP2016101415A (en) * | 2014-11-28 | 2016-06-02 | キヤノン株式会社 | Subject information acquisition apparatus |
| CN105997153A (en) * | 2016-07-15 | 2016-10-12 | 华中科技大学 | Imaging method of ultrasonic CT |
| WO2018116594A1 (en) * | 2016-12-22 | 2018-06-28 | 富士フイルム株式会社 | Photoacoustic image generating device and method, and program |
| CN110772281A (en) * | 2019-10-23 | 2020-02-11 | 哈尔滨工业大学（深圳） | Ultrasonic CT imaging system based on improved ray tracing method |
| WO2020082934A1 (en) * | 2018-10-25 | 2020-04-30 | 南京大学 | Photoacoustic image reconstruction method for inhibiting artifact |
| CN111419185A (en) * | 2020-04-08 | 2020-07-17 | 国网山西省电力公司电力科学研究院 | Magneto-acoustic imaging image reconstruction method with nonuniform sound velocity |
| CN116468859A (en) * | 2023-06-19 | 2023-07-21 | 之江实验室 | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2018061716A (en) * | 2016-10-13 | 2018-04-19 | キヤノン株式会社 | Information processing device, information processing method, and program |
| KR20220069176A (en) * | 2020-11-19 | 2022-05-27 | 삼성전자주식회사 | Photoacoustic apparatus, apparatus and method for obtaining photoacoustic image |
-
2023
- 2023-10-25 CN CN202311387930.2A patent/CN117158911B/en active Active
Patent Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101214156A (en) * | 2008-01-10 | 2008-07-09 | 复旦大学 | Acoustic velocity inhomogeneous medium thermoacoustic imaging reconstruction algorithm |
| JP2012075464A (en) * | 2010-09-30 | 2012-04-19 | Fujifilm Corp | Photoacoustic diagnostic imaging equipment, image generating method, and program |
| CN103445765A (en) * | 2013-09-24 | 2013-12-18 | 南京大学 | Acoustic velocity correction method for photoacoustic imaging |
| WO2015175431A1 (en) * | 2014-05-12 | 2015-11-19 | University Of Washington | Real-time photoacoustic and ultrasound imaging system and method |
| WO2015189268A2 (en) * | 2014-06-10 | 2015-12-17 | Ithera Medical Gmbh | Device and method for hybrid optoacoustic tomography and ultrasonography |
| JP2016101415A (en) * | 2014-11-28 | 2016-06-02 | キヤノン株式会社 | Subject information acquisition apparatus |
| CN104545811A (en) * | 2014-12-26 | 2015-04-29 | 深圳先进技术研究院 | Intravascular imaging system and method |
| CN105249993A (en) * | 2015-11-16 | 2016-01-20 | 南京大学 | Method for selecting optimum sound velocity group to optimize ultrasonic imaging through photoacoustic imaging |
| CN105997153A (en) * | 2016-07-15 | 2016-10-12 | 华中科技大学 | Imaging method of ultrasonic CT |
| WO2018116594A1 (en) * | 2016-12-22 | 2018-06-28 | 富士フイルム株式会社 | Photoacoustic image generating device and method, and program |
| WO2020082934A1 (en) * | 2018-10-25 | 2020-04-30 | 南京大学 | Photoacoustic image reconstruction method for inhibiting artifact |
| CN110772281A (en) * | 2019-10-23 | 2020-02-11 | 哈尔滨工业大学（深圳） | Ultrasonic CT imaging system based on improved ray tracing method |
| CN111419185A (en) * | 2020-04-08 | 2020-07-17 | 国网山西省电力公司电力科学研究院 | Magneto-acoustic imaging image reconstruction method with nonuniform sound velocity |
| CN116468859A (en) * | 2023-06-19 | 2023-07-21 | 之江实验室 | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity |
Non-Patent Citations (4)
| Title |
|---|
| 一种基于光声成像选取最佳声速组的研究;叶濛, 袁杰;《 南京大学学报(自然科学)》;第52卷(第03期);全文 * |
| 光声成像及其在生物医学中的应用;谷怀民;杨思华;向良忠;;生物化学与生物物理进展(05);全文 * |
| 拟合延时补偿的光声重建;朱毅;李文超;张丽君;袁杰;;光电工程(04);全文 * |
| 生物声学成像中声速不均匀性解决方法的研究进展;孙正;贾艺璇;;声学技术(05);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117158911A (en) | 2023-12-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR102144551B1 (en) | Laser optoacoustic ultrasonic imaging system (louis) and methods of use | |
| US10231707B2 (en) | Ultrasound waveform tomography with wave-energy-based preconditioning | |
| US8870770B2 (en) | Low-cost device for C-scan acoustic wave imaging | |
| CN100446730C (en) | Photoacoustic imaging and chromatographic imaging method based on acoustic lens and apparatus thereof | |
| JP5441795B2 (en) | Imaging apparatus and imaging method | |
| JP2015507947A5 (en) | ||
| EP2482713B1 (en) | Photoacoustic measuring apparatus | |
| KR20140103932A (en) | Dual modality imaging system for coregistered functional and anatomical mapping | |
| JPS62148654A (en) | Method and apparatus for enhancing probability of sonic wavemeasurement in living body in organ tissue | |
| CN116468859B (en) | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity | |
| JP2017140092A (en) | Subject information acquisition device | |
| CN117158911B (en) | Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method | |
| Omidvar et al. | Shape estimation of flexible ultrasound arrays using spatial coherence: A preliminary study | |
| Ruiter et al. | Ultrasound Tomography | |
| US20050160817A1 (en) | Superresolution ultrasound | |
| US20240180524A1 (en) | Divergent lens array | |
| CN117320635A (en) | Diverging lens array | |
| Amanatiadis et al. | Transcranial ultrasonic propagation and enhanced brain imaging exploiting the focusing effect of the skull | |
| Qu et al. | Phase aberration correction by bent-ray tracing method for ultrasound computed tomography | |
| JP2018000305A (en) | Subject information acquisition device and signal processing method | |
| Pan et al. | Effect of center frequency and bandwidth of ultrasonic transducer on photoacoustic tomography based on virtual PAT | |
| Zhao et al. | Limited-View Photoacoustic Tomography Based on Transient Radiative Transfer Equation | |
| Niederer | Ultrasound imaging and Doppler flow velocity measurement | |
| Opielinski et al. | Ultrasonic projection imaging of biological media | |
| Thotakura | Simulation of Ultrasound Computed Tomography in Diffraction Mode |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20240401 Address after: Room 1909, Building A, No. 274 Tianmushan Road and No. 2-18 (Shuang) Wantang Road, Xihu District, Hangzhou City, Zhejiang Province, 310012 Patentee after: Zhejiang Liying Medical Technology Co.,Ltd. Country or region after: China Address before: Room 1004, Building 8, Zhihui Zhongchuang Center, Xihu District, Hangzhou City, Zhejiang Province, 310012 Patentee before: Hangzhou Liying Optoelectronic Imaging Co.,Ltd. Country or region before: China |