CN116468859B - Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity - Google Patents
Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity Download PDFInfo
- Publication number
- CN116468859B CN116468859B CN202310722329.8A CN202310722329A CN116468859B CN 116468859 B CN116468859 B CN 116468859B CN 202310722329 A CN202310722329 A CN 202310722329A CN 116468859 B CN116468859 B CN 116468859B
- Authority
- CN
- China
- Prior art keywords
- region
- point
- observation
- sound velocity
- interest
- 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 58
- 238000003325 tomography Methods 0.000 title claims abstract description 56
- 238000009828 non-uniform distribution Methods 0.000 title claims abstract description 30
- 239000002609 medium Substances 0.000 claims description 16
- 238000005070 sampling Methods 0.000 claims description 15
- 238000012216 screening Methods 0.000 claims description 15
- 230000014509 gene expression Effects 0.000 claims description 14
- 239000012736 aqueous medium Substances 0.000 claims description 3
- 239000007787 solid Substances 0.000 claims description 3
- 210000004556 brain Anatomy 0.000 description 17
- 238000004364 calculation method Methods 0.000 description 11
- 238000012545 processing Methods 0.000 description 11
- 239000004973 liquid crystal related substance Substances 0.000 description 10
- 238000010586 diagram Methods 0.000 description 5
- 238000009826 distribution Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000010521 absorption reaction Methods 0.000 description 2
- 238000004590 computer program Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000002604 ultrasonography Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000011664 signaling Effects 0.000 description 1
- 239000012798 spherical particle Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0093—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy
- A61B5/0095—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy by applying light and detecting acoustic waves, i.e. photoacoustic measurements
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Animal Behavior & Ethology (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- Surgery (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Computer Graphics (AREA)
- Acoustics & Sound (AREA)
- Theoretical Computer Science (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
The application discloses a three-dimensional photoacoustic tomography method and a device suitable for non-uniform distribution of sound velocity, wherein the method mainly comprises the following steps: setting corresponding sound velocity areas for an interested area and other areas in the observation view field area respectively; calculating the intersection point of the connecting line of each point to be imaged and each array element in the region of interest in the observation view field region and the region of interest in the observation view field region; calculating the path length of each point to be imaged in the region of interest of the observation field region to reach different sound velocity regions where each array element passes; according to different sound speeds in different sound speed areas and different path lengths, the point (x, y, z) in the observation field area reaches the time delay point required by the ith array element, so that the photoacoustic signal sent by the point (x, y, z) in the observation field area received by the ith array element is calculated, and the photoacoustic signals sent by the point (x, y, z) in the observation field area received by all the array elements are subjected to weighted time delay summation, so that the photoacoustic image reconstructed in the observation field area is obtained.
Description
Technical Field
The application relates to the field of photoacoustic tomography, in particular to a three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity.
Background
Photoacoustic tomography (Photoacoustic computed tomography) has become a widely focused non-destructive medical imaging technique in recent years. The principle of this technique is that in photoacoustic tomography of biological tissue, the tissue is irradiated with electromagnetic waves (usually laser or microwave) pulses and then thermally expanded, thereby exciting photoacoustic waves. Photoacoustic waves are ultrasonic waves that carry the distribution characteristics of absorption of electromagnetic waves by tissue. The electromagnetic wave absorption distribution inside the tissue can be obtained by detecting the photoacoustic wave by the ultrasonic transducer based on a corresponding reconstruction algorithm.
The sound wave propagates in different tissues of an organism, and the propagation sound speed is also different, and in most application researches, the system sound speed is usually set to be constant for imaging. However, the difference in sound speed affects the time of arrival of the photoacoustic signal at the ultrasound transducer, and thus the distance of the sound source from the ultrasound transducer calculated from the propagation time. The assumption of constant sound velocity can lead to errors between the position of the target tissue in the reconstructed image and the actual position of the target tissue, so that problems such as detail blurring, target dislocation and serious artifacts exist in the image.
Therefore, photoacoustic tomography is urgently required for a technique of sound velocity optimization. The currently known sound velocity optimization technique is mainly directed to two-dimensional image reconstruction, even three-dimensional image reconstruction, and is also a pseudo three-dimensional image formed by means of superposition of multiple frames of two-dimensional images. In addition, the sound velocity optimization adopts a complex calculation mode, and the calculation cost is high.
Disclosure of Invention
Aiming at the defects of the prior art, the application provides a three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity. According to a first aspect of an embodiment of the present application, there is provided a three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity, including:
step S1, setting a region of interest in an observation view field region as a first sound velocity region, and setting other regions except the region of interest in the observation view field region as a second sound velocity region; the sound velocity in the first sound velocity zone is the sound velocity corresponding to the tissue to be detected and is recorded as the first sound velocity; the sound velocity of the second sound velocity zone is the sound velocity corresponding to the medium positioned at the periphery of the tissue to be measured, and is recorded as the second sound velocity.
Step S2, calculating the intersection point of each to-be-imaged point and each ultrasonic transducer array element connecting line in the interested region in the observation view field region and the interested region in the observation view field region;
step S3, calculating a first path length of each point to be imaged in the interested region of the observation field region passing through the first sound velocity region and a second path length of each point to be imaged passing through the second sound velocity region;
step S4, calculating a first delay point corresponding to the flight time of the photoacoustic signal in the first sound velocity region in the observation view field region according to the first sound velocity, the first path length and the sampling frequency; calculating a second delay point corresponding to the flight time of the photoacoustic signal in the second sound speed region in the observation view field region according to the second sound speed, the second path length and the sampling frequency; and taking the sum of the first delay point and the second delay point as a delay point required by observing the point in the view field area to the ith array element, calculating the photoacoustic signal emitted by the point in the view field area received by the ith array element, and carrying out weighted delay summation on the photoacoustic signals emitted by the points in the view field area received by all the array elements to obtain a photoacoustic image reconstructed by the view field area.
According to a second aspect of the embodiment of the present application, there is provided a three-dimensional photoacoustic tomography apparatus suitable for non-uniform distribution of sound velocity, including one or more processors for implementing the above-described three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity.
According to a third aspect of the embodiments of the present application, there is provided a computer-readable storage medium having stored thereon a program for implementing the above-described three-dimensional photoacoustic tomography method applicable to non-uniform distribution of sound velocity when executed by a processor.
Compared with the prior art, the application has the beneficial effects that: the application provides a three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity, which reasonably models a non-uniform sound environment in a living body, realizes quick photoacoustic image reconstruction on an observation view field area with non-uniform distribution of sound velocity, effectively solves the problems of detail blurring, target dislocation, serious artifact and the like in a three-dimensional reconstruction image caused by mismatch of reconstruction sound velocity, and can remarkably improve the reconstruction quality of the image. The method can realize three-dimensional photoacoustic image reconstruction with high frame rate and high resolution of the object to be detected under the condition of not increasing hardware overhead and photoacoustic image reconstruction calculation amount. The method provided by the application not only has wide application scenes, but also has more practical application value in the study of photoacoustic tomography technology.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings that are needed in the description of the embodiments will be briefly described below, it being obvious that the drawings in the following description are only some embodiments of the present application, and that other drawings may be obtained according to these drawings without inventive effort to a person skilled in the art.
Fig. 1 is a flowchart of a three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity provided by an embodiment of the present application.
Fig. 2 is a view of the result of XOY plane projection of three-dimensional photoacoustic tomography of mouse brain using constant sound velocity according to an embodiment of the present application.
Fig. 3 is a view of the result of XOY plane projection of three-dimensional photoacoustic tomography of the brain of a mouse using the method of the present application.
Fig. 4 is a view of the result of XOZ plane projection of three-dimensional photoacoustic tomography using constant sound velocity for mouse brain provided by the embodiment of the present application.
Fig. 5 is a view of an XOZ plane projection result of three-dimensional photoacoustic tomography of a mouse brain using the method provided by the embodiment of the present application.
Fig. 6 is a view of the YOZ plane projection result of three-dimensional photoacoustic tomography using constant sound velocity for mouse brain according to the embodiment of the present application.
Fig. 7 is a view of the YOZ plane projection result of three-dimensional photoacoustic tomography of the brain of a mouse according to the embodiment of the present application.
Fig. 8 is a schematic diagram of a three-dimensional photoacoustic tomography method apparatus suitable for non-uniform distribution of sound velocity according to an embodiment of the present application.
Detailed Description
In order that the above objects, features and advantages of the application will be readily understood, a more particular description of the application will be rendered by reference to the appended drawings. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present application. The application may be embodied in many other forms than described herein and similarly modified by those skilled in the art without departing from the spirit or scope of the application, which is therefore not limited to the specific embodiments disclosed below.
The following detailed description of preferred embodiments of the application is made in connection with the accompanying drawings, which form a part hereof, and together with the description of the embodiments of the application, are used to explain the principles of the application and are not intended to limit the scope of the application.
The application aims to provide a three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity, which reasonably models a non-uniform sound environment in a living body while adding no additional calculation cost as much as possible so as to solve the problems of detail blurring, target dislocation, serious artifact and the like in a three-dimensional reconstructed image caused by sound velocity mismatch.
In order to achieve the above purpose, the core of the technical scheme adopted by the application is to set the interested region in the observation view field region as a first sound velocity region according to the position of the interested region in the observation view field, and the rest regions in the observation view field region are second sound velocity regions. And calculating the path lengths of the points in the region of interest to reach different sound velocity regions penetrated by each array element, so as to reconstruct a three-dimensional photoacoustic tomography reconstruction image in the observation field of view based on a delay summation method with a solid angle weighting factor.
As shown in fig. 1, the three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity provided by the embodiment of the application mainly comprises the following steps:
s1, setting a region of interest in an observation view field region as a first sound velocity region, and setting other regions except the region of interest in the observation view field region as a second sound velocity region; the sound velocity in the first sound velocity zone is the sound velocity corresponding to the tissue to be detected and is recorded as the first sound velocity; the sound velocity of the second sound velocity zone is the sound velocity corresponding to the medium positioned at the periphery of the tissue to be measured, and is recorded as the second sound velocity.
And setting an observation field area for reconstructing the three-dimensional photoacoustic tomography image under a Cartesian three-dimensional coordinate system in which the three-dimensional photoacoustic tomography imaging device is located. Wherein the three-dimensional sound velocity region is set so as to be contained therein according to the approximate position of the region of interest in the observation field region. Without loss of generality, considering the shape of tissue in most living bodies, the region of interest in the observation view field region is set to an ellipsoid (which can also be transformed into a sphere depending on the application), and the ellipsoid is a standard ellipsoid.
The medium around the tissue to be measured is a tissue medium or an aqueous medium around the tissue to be measured.
S2, calculating the intersection point of each to-be-imaged point and each ultrasonic transducer array element connecting line in the interested region in the observation view field region and the interested region in the observation view field region.
The step S2 includes: solving the intersection point of a certain point to be imaged in the region of interest of the observation field region and a certain array element connecting line with a preset ellipsoidal region of interest。
Assume that a certain point to be imaged in a region of interest of an observation field area has coordinates ofThe coordinates of certain ultrasonic transducer array element are +.>The spatial straight line through these two points can be expressed as:
assume that the center coordinates of an ellipsoid areThe radius of the shaft is +.>The corresponding standard ellipsoidal equation is:
the solution of the simultaneous straight line equation and the ellipsoid equation is the intersection point of the space straight line and the ellipsoid.
Classification discussions are required when solving equations:
when the x-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are different, screening the intersection points obtained by solving the simultaneous equations according to the x-axis components of the intersection points; the method comprises the following steps:
when (when)In the time-course of which the first and second contact surfaces,
wherein, the liquid crystal display device comprises a liquid crystal display device,for->And (4) point->Slope of the line->Is the intercept of the above-mentioned wire. />Is taken as a pointAnd (4) point->Slope of the line->Is the intercept of the above-mentioned wire.
The following results were obtained by calculation after simultaneous use:
from the following componentsCan calculate the intersection point of the space straight line and the ellipsoid +.>Is +.>The components are as follows:
according to the intersection point and each point to be imaged in the interested area of the observation field areaThe above elements->The calculated coordinates of the intersection point are filtered. The screening can be performed according to the following formula:
it should be noted that since the intersection point is located at each point to be imaged in the region of interest of the observation field regionThe above elements->Between, i.e. intersection->Is +.>The component size is located in the above viewEach point to be imaged in the region of interest of the region of view>Is +.>Component and the above-mentioned array element->Is +.>Between the components, the intersection point can be screened out by the two formulas>Is +.>The effective value of the component.
From the crossing points after screeningIs +.>The components can be calculated to obtain the intersection +.>Is +.>Component sum->A component.
When the x-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are the same and the y-axis coordinates are different, screening intersection points obtained by simultaneous equation solving according to y-axis components of the intersection points; the method comprises the following steps:
when (when)Time, emptyThe straight line between->On the plane of (a), the problem is converted into solutionThe intersection point of the upper plane straight line and the ellipsoid.
Wherein, the liquid crystal display device comprises a liquid crystal display device,for->And (4) point->Slope of the line->Is the intercept of the above-mentioned wire.
The following results were obtained by calculation after simultaneous use:
from the following componentsCan calculate the intersection point of the space straight line and the ellipsoid +.>Is +.>The components are as follows:
based on the intersection point and the observation field regionEach point to be imaged in the interest areaThe above elements->The calculated coordinates of the intersection point are filtered. The screening can be performed according to the following formula:
it should be noted that since the intersection point is located at each point to be imaged in the region of interest of the observation field regionThe above elements->Between, i.e. intersection->Is +.>The component size is within the region of interest of the above-mentioned observation field region and the point to be imaged is +.>Is +.>Component and the above-mentioned array element->Is +.>Between the components, the intersection point can be screened out by the two formulas>Is +.>The effective value of the component.
From the crossing points after screeningIs +.>The components can be calculated to obtain the intersection +.>Is +.>A component.
And when the x-axis coordinates and the y-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are the same and the z-axis coordinates are different, screening the intersection points obtained by simultaneous solving according to the z-axis components of the intersection points. The method comprises the following steps:
when (when)The spatial line is a line parallel to the z-axis and is co-ordinate +.>As is known, the following equation is followed:
from the following componentsCan calculate the intersection point of the space straight line and the ellipsoid +.>Is +.>The components are as follows:
according to the intersection point and the point to be imaged in the interested area of the observation field areaThe above elements->The calculated coordinates of the intersection point are filtered. The screening can be performed according to the following formula:
。
it should be noted that since the intersection point is located at each point to be imaged in the region of interest of the observation field regionThe above elements->Between, i.e. intersection->Is +.>The component size is within the region of interest of the above-mentioned observation field region and the point to be imaged is +.>Is +.>Component and the above-mentioned array element->Is +.>Between the components, the intersection point can be screened out by the two formulas>Is +.>The effective value of the component.
S3, calculating a first path length of each point to be imaged in the region of interest of the observation field region, which passes through the first sound velocity region, and a second path length of each point to be imaged, which passes through the second sound velocity region.
From the previous calculation results, the point to be imaged in the interested area of the observation field area can be calculatedCrossing point->The distance between the two is obtained to obtain a first path length passing through a first sound velocity zone, and the expression is as follows:
and the point to be imaged in the region of interestIs combined with the array element->The distance between them, giving the total length, is expressed as follows:
the array elementCrossing point->Distance between->It is readily available to calculate the difference between the total length and the first path length to obtain a second path length through the second sound velocity zone as follows:
in order to improve the reconstruction effect of the region of interest in the observation field region and to screen out the complex solutions introduced in the remaining regions of the observation field region in the above calculation process, the calculated distance distribution values in the remaining regions of the observation field region need to be screened and processed according to the following equation set.
。
The meaning of the two formulas above is: at a distance ofAnd distance->The points satisfying the formula in brackets are located outside the ellipsoidal region of interest, since they are +.>The connecting line and the ellipsoidal interested area have no intersection point, so that complex solutions introduced in the rest of the observation field area in the calculation process can be screened out by eliminating the connecting line and the ellipsoidal interested area.
S4, calculating a first delay point corresponding to the flight time of the photoacoustic signal in the first sound velocity region in the observation view field region according to the first sound velocity, the first path length and the sampling frequency; calculating a second delay point corresponding to the flight time of the photoacoustic signal in the second sound speed region in the observation view field region according to the second sound speed, the second path length and the sampling frequency; taking the sum of the first time delay point and the second time delay point as a point in the observation view field areaTo->The delay point required by each array element is calculated by this>Point in the observation view field region received by the individual array elements>The emitted photoacoustic signal is about the point in the observation view field area received by all array elements>And carrying out weighted delay summation on the sent photoacoustic signals to obtain a photoacoustic image reconstructed in the observation field area.
The expression of the photoacoustic image reconstructed by calculating the observation field region is as follows:
wherein, the liquid crystal display device comprises a liquid crystal display device,indicating the point in the region of the field of view +.>Is used for the image intensity values of (a),indicating->The array elements are directed to the point in the observation field area>The weighting factor corresponding to the opening angle is used,indicate->The observation field area points received by each array element>The emitted photo-acoustic signal is used to generate,indicating the point in the region of the field of view +.>To->The delay points required by the array elements can be expressed as:
wherein, the liquid crystal display device comprises a liquid crystal display device,and the time delay points corresponding to the flight time of the photoacoustic signals in different areas of the indicated observation field.
Let the sound velocity of the region of interest in the observation field region beThe sound velocity of the rest of the observation view field region is。
Then there are:
wherein, the liquid crystal display device comprises a liquid crystal display device,represents the sampling frequency of an ultrasonic transducer array element, +.>Represented is a rounding down.
Example 1
The embodiment of the application describes three-dimensional photoacoustic image reconstruction of the brain of a mouse by using a bowl-shaped array through the method. Wherein 1024 array elements are uniformly distributed on the inner surface of the bowl-shaped array, the radius of the bowl-shaped array is 100mm, and the sampling rate is high。
Referring next to fig. 1, photoacoustic image reconstruction is performed using the method provided by the present patent according to the following steps:
s1, setting a region of interest in an observation view field region as a first sound velocity region, and setting other regions except the region of interest in the observation view field region as a second sound velocity region; the sound velocity in the first sound velocity zone is the sound velocity corresponding to the tissue to be detected and is recorded as the first sound velocity; the sound velocity of the second sound velocity zone is the sound velocity corresponding to the medium positioned at the periphery of the tissue to be measured, and is recorded as the second sound velocity.
Under a Cartesian three-dimensional coordinate system where the three-dimensional photoacoustic tomography experiment device is located, according to the position of the head of a mouse in the experiment, a cube with observation view field areas of 20mm multiplied by 20mm for reconstructing the three-dimensional photoacoustic tomography image is set, and the imaging resolution is 50 mu m. The region of interest in the observation view field region, i.e., the region where the mouse brain is located, is set to be a standard ellipsoid, expressed as:
wherein the center coordinate of the ellipsoid isThe axial radii of the ellipsoids are respectively. The sound velocity in this region is the corresponding sound velocity of the mouse brain. The remaining region sound speed is the corresponding sound speed of the aqueous medium.
S2, calculating the intersection point of each to-be-imaged point and each ultrasonic transducer array element connecting line in the interested region in the observation view field region and the interested region in the observation view field region.
The intersection point of a certain point to be imaged in the region of interest of the observation field area and a certain array element connecting line with a preset ellipsoidal region of interest is set as. Selecting single coordinates in a matrix of points to be imaged in a region of interest of an observation field of view one by one with a double loop>And a single coordinate in the matrix of array elements +.>The spatial straight line through these two points can be expressed as:
and (2) combining the linear equation and the ellipsoidal equation established in the step (S1), wherein the solution is the intersection point of the space linear and the ellipsoid.
Classification discussions are required when solving equations:
when the x-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are different, screening the intersection points obtained by solving the simultaneous equations according to the x-axis components of the intersection points;
when the x-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are the same and the y-axis coordinates are different, screening intersection points obtained by simultaneous equation solving according to y-axis components of the intersection points;
and when the x-axis coordinates and the y-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are the same and the z-axis coordinates are different, screening the intersection points obtained by solving the simultaneous equations according to the z-axis components of the intersection points.
S3, calculating a first path length of each point to be imaged in the region of interest of the observation field region, which passes through the first sound velocity region and a second path length of each point to be imaged, which passes through the second sound velocity region.
Calculating the distance between a certain point to be imaged in the region of interest of the observation field region and the intersection point to obtain a first path length passing through a first sound velocity region;
calculating the distance between a certain point to be imaged and a certain ultrasonic transducer array element in an interested region of the observation field region to obtain the total length;
and calculating the difference between the total length and the first path length to obtain a second path length passing through the second sound speed zone.
Further, the step S3 further includes screening the first path length and the second path length, where the expression is as follows:
；
in the method, in the process of the application,for the first path length, +>For the second path length->For observing the coordinates of a certain point to be imaged in the region of interest of the field of view>To observe the ellipsoidal central coordinates of a region of interest in the field of view,the axial radii of the ellipsoids, respectively.
S4, calculating a first delay point corresponding to the flight time of the photoacoustic signal in the first sound velocity region in the observation view field region according to the first sound velocity, the first path length and the sampling frequency; calculating a second delay point corresponding to the flight time of the photoacoustic signal in the second sound speed region in the observation view field region according to the second sound speed, the second path length and the sampling frequency; and taking the sum of the first delay point and the second delay point as a delay point required by observing the point in the view field area to the ith array element, calculating the photoacoustic signal emitted by the point in the view field area received by the ith array element, and carrying out weighted delay summation on the photoacoustic signals emitted by the points in the view field area received by all the array elements to obtain a photoacoustic image reconstructed by the view field area.
Specifically, the step S4 includes:
calculating a first delay point corresponding to the flight time of the photoacoustic signal in the first sound velocity region in the observation view field region according to the first sound velocity, the first path length and the sampling frequency of the ultrasonic transducer array element;
calculating a second delay point corresponding to the flight time of the photoacoustic signal in the second sound speed region in the observation view field region according to the second sound speed, the second path length and the sampling frequency of the ultrasonic transducer array element;
taking the sum of the first time delay point and the second time delay point as a point in the observation view field areaTo the delay point required by the ith array element;
from points within the field of viewCalculating the time delay point required from the ith array element to the ith array elementThe observation field of view area internal point received by each array elementAn emitted photoacoustic signal;
setting the firstThe array elements are directed to the points in the observation field areaA weighting factor corresponding to the opening angle;
delay and sum method based on weighting factors with solid anglesPoints in the observation field of view area received for all array elementsThe sent photoacoustic signals are combined with the weighting factors corresponding to the opening angles to carry out weighted summation, so that the point +.>And obtaining a photoacoustic image reconstructed from the observation field region.
Fig. 2 shows an XOY plane projection result diagram of three-dimensional photoacoustic tomography with constant sound velocity for the mouse brain. Fig. 3 shows an XOY plane projection result diagram of three-dimensional photoacoustic tomography of the brain of a mouse in example 1. Fig. 4 shows an XOZ plane projection result graph of three-dimensional photoacoustic tomography with constant sound velocity for the brain of a mouse. Fig. 5 shows an XOZ plane projection result diagram of three-dimensional photoacoustic tomography of the brain of a mouse of example 1. Fig. 6 shows a view of the YOZ plane projection result of three-dimensional photoacoustic tomography with constant sound velocity for the brain of a mouse. Fig. 7 shows a YOZ plane projection result graph of three-dimensional photoacoustic tomography of the mouse brain using the method provided in example 1.
Referring to fig. 2 to 7, by comparing each plane projection result of three-dimensional photoacoustic tomography by directly using a constant sound velocity on a mouse brain with each plane projection result of three-dimensional photoacoustic tomography by using an embodiment, it can be seen that, since there is a significant difference between the sound velocity in a mouse brain region and a surrounding water medium region, each plane projection result of three-dimensional photoacoustic tomography by directly using a constant sound velocity is unclear, and the reconstructed image has problems of characteristic divergence, detail ambiguity, target dislocation, artifact, and the like. This is caused by the mismatch of sound speeds during reconstruction. The projection results of each plane of the three-dimensional photoacoustic tomography by adopting the method provided by the application have higher definition, and the detail characteristics of the reconstruction region can be accurately reflected. The embodiment of the application is accurate in compensating the sound velocity distribution in the observation view field area, and the method provided by the application has a good practical application effect.
In summary, the three-dimensional photoacoustic tomography method suitable for the non-uniform distribution of the sound velocity provided by the application can reasonably model the non-uniform acoustic environment in the living body under the condition of not additionally increasing the system preprocessing and calculation cost, realize the rapid photoacoustic image reconstruction on the observation view field area with the non-uniform distribution of the sound velocity, effectively solve the problems of detail blurring, target dislocation, serious artifact and the like in the three-dimensional reconstruction image caused by the mismatch of the reconstructed sound velocity, and can remarkably improve the reconstruction quality of the image. The method can realize three-dimensional photoacoustic image reconstruction with high frame rate and high resolution of the object to be detected under the condition of not increasing hardware overhead and photoacoustic image reconstruction calculation amount. The method provided by the application not only has wide application scenes, but also has more practical application value in the study of photoacoustic tomography technology.
Referring to fig. 8, a three-dimensional photoacoustic tomography method apparatus suitable for non-uniform distribution of sound velocity according to an embodiment of the present application includes one or more processors for implementing the three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity in the above embodiment.
The embodiment of the application based on the three-dimensional photoacoustic tomography method device suitable for the non-uniform distribution of sound velocity can be applied to any device with data processing capability, and the device with data processing capability can be a device or a device such as a computer. The apparatus embodiments may be implemented by software, or may be implemented by hardware or a combination of hardware and software. Taking software implementation as an example, the device in a logic sense is formed by reading corresponding computer program instructions in a nonvolatile memory into a memory by a processor of any device with data processing capability. In terms of hardware, as shown in fig. 8, a hardware structure diagram of an apparatus with optional data processing capability where the photoacoustic signaling device based on the spherical particle light pulse excitation effect of the present application is located is shown in fig. 8, and in addition to the processor, the memory, the network interface, and the nonvolatile memory shown in fig. 8, the apparatus with optional data processing capability where the apparatus with optional data processing capability is located in the embodiment generally includes other hardware according to the actual function of the apparatus with optional data processing capability, which is not described herein again.
The implementation process of the functions and roles of each unit in the above device is specifically shown in the implementation process of the corresponding steps in the above method, and will not be described herein again.
For the device embodiments, reference is made to the description of the method embodiments for the relevant points, since they essentially correspond to the method embodiments. The apparatus embodiments described above are merely illustrative, wherein the elements illustrated as separate elements may or may not be physically separate, and the elements shown as elements may or may not be physical elements, may be located in one place, or may be distributed over a plurality of network elements. Some or all of the modules may be selected according to actual needs to achieve the purposes of the present application. Those of ordinary skill in the art will understand and implement the present application without undue burden.
The embodiment of the present application also provides a computer-readable storage medium having stored thereon a program which, when executed by a processor, implements the three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity in the above-described embodiment.
The computer readable storage medium may be an internal storage unit, such as a hard disk or a memory, of any of the data processing enabled devices described in any of the previous embodiments. The computer readable storage medium may be any device having data processing capability, for example, a plug-in hard disk, a Smart Media Card (SMC), an SD Card, a Flash memory Card (Flash Card), or the like, which are provided on the device. Further, the computer readable storage medium may include both internal storage units and external storage devices of any data processing device. The computer readable storage medium is used for storing the computer program and other programs and data required by the arbitrary data processing apparatus, and may also be used for temporarily storing data that has been output or is to be output.
It will be appreciated by persons skilled in the art that the foregoing description is a preferred embodiment of the application, and is not intended to limit the application, but rather to limit the application to the specific embodiments described, and that modifications may be made to the technical solutions described in the foregoing embodiments, or equivalents may be substituted for elements thereof, for the purposes of those skilled in the art. Modifications, equivalents, and alternatives falling within the spirit and principles of the application are intended to be included within the scope of the application.
Claims (10)
1. A three-dimensional photoacoustic tomography method suitable for non-uniform distribution of sound velocity, comprising:
step S1, setting a region of interest in an observation view field region as a first sound velocity region, and setting other regions except the region of interest in the observation view field region as a second sound velocity region; the sound velocity in the first sound velocity zone is the sound velocity corresponding to the tissue to be detected and is recorded as the first sound velocity; the sound velocity of the second sound velocity zone is the sound velocity corresponding to the medium positioned at the periphery of the tissue to be detected and is recorded as the second sound velocity;
step S2, calculating the intersection point of each to-be-imaged point and each ultrasonic transducer array element connecting line in the interested region in the observation view field region and the interested region in the observation view field region;
step S3, calculating a first path length of each point to be imaged in the interested region of the observation field region passing through the first sound velocity region and a second path length of each point to be imaged passing through the second sound velocity region;
step S4, calculating a first delay point corresponding to the flight time of the photoacoustic signal in the first sound velocity region in the observation view field region according to the first sound velocity, the first path length and the sampling frequency; calculating a second delay point corresponding to the flight time of the photoacoustic signal in the second sound speed region in the observation view field region according to the second sound speed, the second path length and the sampling frequency; taking the sum of the first delay point and the second delay point as a delay point required by observing the point in the view field area to the ith array element, calculating a photoacoustic signal emitted by the point in the view field area received by the ith array element, and carrying out weighted delay summation on photoacoustic signals emitted by the points in the view field area received by all the array elements to obtain a photoacoustic image reconstructed by the view field area;
the expressions of the first time delay point and the second time delay point are as follows:
wherein n is 1 (x, y, z, i) is the first delay point, n 2 (x, y, z, i) is the second delay point, d 1 (x, y, z, i) is the first path length, d 2 (x, y, z, i) is the second path length, c 1 At a first sound velocity, c 2 Is the second sound velocity, f s Representing the sampling frequency of the ultrasonic transducer elements,represented is a downward rounding;
wherein n (x, y, z, i) represents the time delay point required for observing the point (x, y, z) from the ith array element in the field of view, n k And (x, y, z, i) representing time delay points corresponding to the flight time of the photoacoustic signals in different areas of the observation field of view.
2. A three-dimensional photoacoustic tomography method suitable for use in a non-uniform distribution of sound speeds according to claim 1, wherein the region of interest in the region of view field is ellipsoidal or spherical.
3. The three-dimensional photoacoustic tomography method of claim 1, wherein the medium around the tissue to be measured is a tissue medium or an aqueous medium around the tissue to be measured.
4. A three-dimensional photoacoustic tomography method for non-uniform distribution of sound speed according to claim 1, wherein calculating the intersection point of the connection line of each image point to be imaged and each ultrasonic transducer array element in the region of interest in the observation view field region and the region of interest in the observation view field region comprises:
acquiring coordinates of a certain to-be-imaged point and coordinates of a certain ultrasonic transducer array element in an interested region of an observation field region, and acquiring an expression of a connecting line of the certain to-be-imaged point and the certain ultrasonic transducer array element in the interested region of the observation field region;
acquiring the ellipsoidal central coordinate and the axial radius of the region of interest in the observation view field region as well as the expression of the region of interest in the observation view field region;
and solving an expression of a connecting line of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the simultaneous observation view field region and an expression of the region of interest of the observation view field region to obtain an intersection point.
5. The method of claim 4, wherein the step of solving the intersection point by using the expression of the line connecting the point to be imaged and the ultrasonic transducer element in the region of interest of the view field and the expression of the region of interest of the view field, together, comprises the steps of:
when the x-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are different, screening the intersection points obtained by solving the simultaneous expression according to the x-axis components of the intersection points;
when the x-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are the same and the y-axis coordinates are different, screening intersection points obtained by solving the simultaneous expression according to y-axis components of the intersection points;
when the x-axis coordinates and the y-axis coordinates of a certain to-be-imaged point and a certain ultrasonic transducer array element in the region of interest of the observation field area are the same and the z-axis coordinates are different, the intersection points obtained by solving the simultaneous expression are screened according to the z-axis components of the intersection points.
6. A three-dimensional photoacoustic tomography method suitable for use in a non-uniform distribution of sound speeds according to claim 1 or 4, wherein the step S3 comprises:
calculating the distance between a certain point to be imaged in the region of interest of the observation field region and the intersection point to obtain a first path length passing through a first sound velocity region;
calculating the distance between a certain point to be imaged and a certain ultrasonic transducer array element in an interested region of the observation field region to obtain the total length;
and calculating the difference between the total length and the first path length to obtain a second path length passing through the second sound speed zone.
7. A three-dimensional photoacoustic tomography method for non-uniform distribution of sound speeds according to claim 6, wherein said step S3 further comprises:
screening the first path length and the second path length, wherein the expression is as follows:
wherein d 1 For a first path length, d 2 For a second path length, (x) 1 ,y 1 ,z 1 ) To observe the coordinates of a point to be imaged in the view field region (x) 0 ,y 0 ,z 0 ) For viewing the center coordinates of the ellipsoids of the region of interest in the view field region, a, b, c are the axial radii of the ellipsoids, respectively.
8. A three-dimensional photoacoustic tomography method for non-uniform distribution of sound velocity according to claim 1, wherein the step S4 comprises:
calculating a first delay point corresponding to the flight time of the photoacoustic signal in the first sound velocity region in the observation view field region according to the first sound velocity, the first path length and the sampling frequency of the ultrasonic transducer array element;
calculating a second delay point corresponding to the flight time of the photoacoustic signal in the second sound speed region in the observation view field region according to the second sound speed, the second path length and the sampling frequency of the ultrasonic transducer array element;
taking the sum of the first delay point and the second delay point as a delay point required for observing points (x, y, z) in the view field area to the ith array element;
calculating photoacoustic signals sent by points (x, y, z) in the observation field area received by the ith array element according to time delay points required from the points (x, y, z) in the observation field area to the ith array element;
setting a weighting factor corresponding to the opening angle of the ith array element to the point (x, y, z) in the observation field area;
and carrying out weighted summation on photoacoustic signals sent by points (x, y, z) in the observation field of view area and received by all array elements based on a delay summation method with solid angle weighting factors, and combining the weighting factors corresponding to the opening angles to obtain image intensity values of the points (x, y, z) in the observation field of view area, so as to obtain a photoacoustic image reconstructed in the observation field of view area.
9. A three-dimensional photoacoustic tomography apparatus adapted for non-uniform distribution of sound velocity comprising one or more processors for implementing the three-dimensional photoacoustic tomography method adapted for non-uniform distribution of sound velocity according to any one of claims 1-8.
10. A computer-readable storage medium having stored thereon a program for implementing the three-dimensional photoacoustic tomography method for non-uniform distribution of sound speeds according to any one of claims 1 to 8 when executed by a processor.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310722329.8A CN116468859B (en) | 2023-06-19 | 2023-06-19 | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310722329.8A CN116468859B (en) | 2023-06-19 | 2023-06-19 | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116468859A CN116468859A (en) | 2023-07-21 |
| CN116468859B true CN116468859B (en) | 2023-09-08 |
Family
ID=87179249
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310722329.8A Active CN116468859B (en) | 2023-06-19 | 2023-06-19 | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116468859B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117158911B (en) * | 2023-10-25 | 2024-01-23 | 杭州励影光电成像有限责任公司 | Multi-sound-velocity self-adaptive photoacoustic tomography image reconstruction method |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5384912A (en) * | 1987-10-30 | 1995-01-24 | New Microtime Inc. | Real time video image processing system |
| CN104203110A (en) * | 2012-03-26 | 2014-12-10 | 毛伊图像公司 | Systems and methods for improving ultrasound image quality by applying weighting factors |
| CN105249993A (en) * | 2015-11-16 | 2016-01-20 | 南京大学 | Method for selecting optimum sound velocity group to optimize ultrasonic imaging through photoacoustic imaging |
| CN108364326A (en) * | 2018-02-08 | 2018-08-03 | 中南大学 | A kind of CT imaging methods |
| CN111214213A (en) * | 2020-02-13 | 2020-06-02 | 南京科技职业学院 | Photoacoustic tomography method suitable for medium with nonuniform sound velocity |
| CN113936069A (en) * | 2021-09-29 | 2022-01-14 | 之江实验室 | Array element virtual interpolation method for photoacoustic tomography |
| CN115177217A (en) * | 2022-09-09 | 2022-10-14 | 之江实验室 | Photoacoustic signal simulation method and device based on spherical particle light pulse excitation effect |
| CN115363528A (en) * | 2022-08-03 | 2022-11-22 | 中国科学院深圳先进技术研究院 | Reconstruction method for elevation angle resolution of annular array photoacoustic computed tomography system |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8291348B2 (en) * | 2008-12-31 | 2012-10-16 | Hewlett-Packard Development Company, L.P. | Computing device and method for selecting display regions responsive to non-discrete directional input actions and intelligent content analysis |
| JP7013215B2 (en) * | 2017-11-24 | 2022-01-31 | キヤノン株式会社 | Information processing equipment, information processing method, program |
-
2023
- 2023-06-19 CN CN202310722329.8A patent/CN116468859B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5384912A (en) * | 1987-10-30 | 1995-01-24 | New Microtime Inc. | Real time video image processing system |
| CN104203110A (en) * | 2012-03-26 | 2014-12-10 | 毛伊图像公司 | Systems and methods for improving ultrasound image quality by applying weighting factors |
| CN105249993A (en) * | 2015-11-16 | 2016-01-20 | 南京大学 | Method for selecting optimum sound velocity group to optimize ultrasonic imaging through photoacoustic imaging |
| CN108364326A (en) * | 2018-02-08 | 2018-08-03 | 中南大学 | A kind of CT imaging methods |
| CN111214213A (en) * | 2020-02-13 | 2020-06-02 | 南京科技职业学院 | Photoacoustic tomography method suitable for medium with nonuniform sound velocity |
| CN113936069A (en) * | 2021-09-29 | 2022-01-14 | 之江实验室 | Array element virtual interpolation method for photoacoustic tomography |
| CN115363528A (en) * | 2022-08-03 | 2022-11-22 | 中国科学院深圳先进技术研究院 | Reconstruction method for elevation angle resolution of annular array photoacoustic computed tomography system |
| CN115177217A (en) * | 2022-09-09 | 2022-10-14 | 之江实验室 | Photoacoustic signal simulation method and device based on spherical particle light pulse excitation effect |
Non-Patent Citations (1)
| Title |
|---|
| A Review of High-Frequency Ultrasonic Transducers for Photoacoustic Imaging Applications;Danyang Ren;《IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ( Volume: 69, Issue: 6, June 2022)》;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116468859A (en) | 2023-07-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US20210321874A1 (en) | Transcranial photoacoustic/thermoacoustic tomography brain imaging informed by adjunct image data | |
| US20180228471A1 (en) | Method and apparatus for analyzing elastography of tissue using ultrasound waves | |
| Ding et al. | Model-based reconstruction of large three-dimensional optoacoustic datasets | |
| US9360551B2 (en) | Object information acquiring apparatus and control method thereof | |
| US10531798B2 (en) | Photoacoustic information acquiring apparatus and processing method | |
| CN115177217B (en) | Photoacoustic signal simulation method and device based on spherical particle light pulse excitation effect | |
| US9572531B2 (en) | Object information acquiring apparatus and control method thereof | |
| CN107708575A (en) | Method, system and computer program product for single tracing positional shearing wave elastogram | |
| CN116468859B (en) | Three-dimensional photoacoustic tomography method and device suitable for non-uniform distribution of sound velocity | |
| EP2954839A1 (en) | Photoacoustic apparatus, signal processing method of a photoacoustic apparatus, and program | |
| Rohling | 3D freehand ultrasound: reconstruction and spatial compounding | |
| JP5847454B2 (en) | Subject information acquisition apparatus, display control method, and program | |
| Singh et al. | Synthetic models of ultrasound image formation for speckle noise simulation and analysis | |
| Hakakzadeh et al. | A spatial-domain factor for sparse-sampling circular-view photoacoustic tomography | |
| Duan et al. | Validation of optical-flow for quantification of myocardial deformations on simulated RT3D ultrasound | |
| JP6072206B2 (en) | SUBJECT INFORMATION ACQUISITION DEVICE AND DISPLAY METHOD | |
| JP7275261B2 (en) | 3D ULTRASOUND IMAGE GENERATING APPARATUS, METHOD, AND PROGRAM | |
| Insana et al. | Combining spatial registration with clutter filtering for power-Doppler imaging in peripheral perfusion applications | |
| WO2015174082A1 (en) | Photoacoustic apparatus | |
| Ruiter et al. | P3a-2 resolution assessment of a 3d ultrasound computer tomograph using ellipsoidal backprojection | |
| US20240115234A1 (en) | Ultrasound measurement system and method | |
| JP2018192309A (en) | Information processing device, information processing method, and program | |
| WO2023047601A1 (en) | Image generation method, image generation program, and image generation apparatus | |
| Amanatiadis et al. | Transcranial ultrasonic propagation and enhanced brain imaging exploiting the focusing effect of the skull | |
| JP7479684B2 (en) | IMAGE GENERATION METHOD, IMAGE GENERATION PROGRAM, AND IMAGE GENERATION DEVICE |
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 |