CN109091109A - The image reconstructing method of Optimization-type Photoacoustic tomography based on complete matrix filtering and time reversal operator - Google Patents
The image reconstructing method of Optimization-type Photoacoustic tomography based on complete matrix filtering and time reversal operator Download PDFInfo
- Publication number
- CN109091109A CN109091109A CN201810706816.4A CN201810706816A CN109091109A CN 109091109 A CN109091109 A CN 109091109A CN 201810706816 A CN201810706816 A CN 201810706816A CN 109091109 A CN109091109 A CN 109091109A
- Authority
- CN
- China
- Prior art keywords
- matrix
- information
- filtering
- wave
- imaging
- 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.)
- Granted
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 181
- 238000001914 filtration Methods 0.000 title claims abstract description 49
- 238000000034 method Methods 0.000 title claims abstract description 37
- 238000003325 tomography Methods 0.000 title claims abstract description 20
- 238000002604 ultrasonography Methods 0.000 claims abstract description 17
- 239000006096 absorbing agent Substances 0.000 claims abstract description 6
- 238000003384 imaging method Methods 0.000 claims description 67
- 230000008569 process Effects 0.000 claims description 11
- 230000005540 biological transmission Effects 0.000 claims description 6
- 238000010276 construction Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 5
- 230000017105 transposition Effects 0.000 claims description 4
- 238000010895 photoacoustic effect Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 2
- 238000001514 detection method Methods 0.000 abstract description 2
- 230000009467 reduction Effects 0.000 abstract description 2
- 229910000831 Steel Inorganic materials 0.000 description 4
- 230000008901 benefit Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 239000010959 steel Substances 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 239000004615 ingredient Substances 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 210000004027 cell Anatomy 0.000 description 1
- 239000010432 diamond Substances 0.000 description 1
- 238000011503 in vivo imaging Methods 0.000 description 1
- 230000031700 light absorption Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 210000003463 organelle Anatomy 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000007493 shaping process Methods 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
- 239000011800 void material Substances 0.000 description 1
Classifications
-
- 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
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Signal Processing (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Power Engineering (AREA)
- Acoustics & Sound (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physiology (AREA)
- Psychiatry (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
- Image Processing (AREA)
Abstract
The invention discloses a kind of image reconstructing methods of Photoacoustic tomography based on complete matrix filtering and time reversal operator, first, information matrix is constructed according to the photoacoustic signal that ultrasound transducer array receives, not only included the information containing direct wave in information matrix, but also has included the information of the scattered wave interfered；Then, the inverse diagonal line property according to specific to direct wave part can filter out the scattering wave component in information matrix by operations such as 45 ° of 45 ° of matrix rotations, matrix filtering, matrix counter-rotating reduction；Finally, the spatial distribution image of the emitted intensity of acoustic wave of absorber of light can be obtained from filtered information matrix using time reversal operator.Method proposed by the present invention can effective filter out the interference of scattered wave, while completely remain effective detection information, to effectively improve in scattering medium and photoacoustic image reconstruction accuracy, ease for use and wide applicability with higher after scattering medium.
Description
Technical field
The present invention relates to a kind of images of Optimization-type Photoacoustic tomography based on complete matrix filtering and time reversal operator
Reconstructing method belongs to Photoacoustic tomography technology.
Background technique
A kind of biomedical imaging technology of the Photoacoustic tomography as Noninvasive, had obtained hair by leaps and bounds in recent years
Exhibition.Photoacoustic tomography had not only had the advantages that acoustic method to deep tissues high resolution, but have the function of optical means at
The advantage [1] of picture, molecular imaging and image contrast etc..
However, Photoacoustic tomography is usually applicable in the uniform medium of acoustic characteristic, the imaging of scattering medium is still optoacoustic
The a great problem that imaging technique faces.In scattering medium, in the photoacoustic signal that receives, in addition to comprising from imageable target
Direct wave, also comprising via the scattering wave component after scattering layer Multiple Scattering.Due to the irregular distribution of scatterer, scattered wave
Randomness is presented in ingredient, it is difficult to extract image information, scattered wave can bring speckle noise, pattern distortion to photoacoustic image.When scattered
When ejected wave ingredient occupies leading position, strong scattering wave interference will seriously reduce photoacoustic imaging depth and image quality, even make it
Cognizable image cannot be accessed.
In order to solve this problem, there has been proposed the various photoacoustic imaging schemes suitable for scattering medium.For example, utilizing
It has been proposed that using correlation matrix filtering method promoted scattering medium in photoacoustic imaging quality [2], however this method due to
The photoacoustic data detected cannot be made full use of, the limitations such as that there are imaging regions is narrow, resolution ratio is low.
Bibliography
【1】L.V.Wang and S.Hu,“Photoacoustic tomography:in vivo imaging from
organelles to organs,”Science 335(6075),1458–1462(2012)
【2】Wei Rui,Chao Tao,Xiaojun Liu,“Photoacoustic imaging in scattering
media by combining a correlation matrix filter with a time reversal
operator,”Optics Express 25(19),22840.
Summary of the invention
Goal of the invention: the picture quality in order to promote Photoacoustic tomography in sound scattering medium, the present invention provide a kind of base
In the image reconstructing method of matrix filtering and the Optimization-type Photoacoustic tomography of time reversal operator, effectively to promote sound scattering
The picture quality of photoacoustic imaging in medium.The inverse diagonal line characteristic of information matrix specifically according to corresponding to direct sound wave,
The operation such as reduction is rotated clockwise by 45 ° of matrix rotations counterclockwise, zero padding, matrix filtering, 45 ° of matrix, filters out scattered wave
Then interference signal carries out time reversal imaging, to reduce the interference of scattered sound waves, promote the photoacoustic imaging in scattering medium
Picture quality.
Technical solution: to achieve the above object, the technical solution adopted by the present invention are as follows:
The image reconstructing method of Photoacoustic tomography based on complete matrix filtering and time reversal operator, including walk as follows
It is rapid:
Step 1: using pulsed laser irradiation by Imaged samples, the absorber of light in sample absorbs the energy of pulse laser
Afterwards, ultrasonic wave is radiated to surrounding media due to optoacoustic effect；
Step 2: receiving the ultrasonic signal that sample emission goes out, the ultrasound transducer array using ultrasound transducer array
Transducer unit with the arrangement of N number of straight line along the x-axis direction, the center spacing of two neighboring transducer unit are w, are changed for n-th
The lateral coordinates of energy device unit are as follows:
Step 3: the ultrasonic signal received is carried out including operation time gated, including Short Time Fourier Transform,
Construct the information matrix K that size is N × N；
Step 4: to information matrix K sequence carry out include 45 ° of rotations counterclockwise, zero padding, filtering, 45 ° rotate clockwise and also
Operation including original filters out scattering wave component, retains through wave component；
Step 5: singular value decomposition being carried out to filtered information matrix, to feature corresponding to maximum M singular value
Vector application time inverts operator, obtains the photoacoustic image of corresponding depth；
Step 6: according to the acoustic transfer function in homogeneous medium, acquiring corrected parameter, reduce since limited perspective ultrasound is believed
Number detection caused by including pattern distortion, falseness contrast including image deflects；
Step 7: changing time gated region, repeat step 3~6, obtain the photoacoustic imaging of different depth, finally obtain whole
The photoacoustic image in a region.
Specifically, the construction process of information matrix K is as follows in the step 3:
Step 31: what interception ultrasonic transducer received is inSignal in time window, and it is right
It makees Fourier transformation, obtains the frequency domain presentation P of n-th of transducer unitn(T, f), frequency domain presentation Pn(T, f) includes through
Wave signal PD(T, f) and scattered wave signal PSTwo parts (T, f)；Wherein, T indicates that imaging plane is passed to transducer array sound wave
It is the shortest time broadcast, also referred to as time gated；Δ t indicates that the width of time window, f indicate frequency；
So-called direct wave refers to from imaging region, without the refraction of any scatterer, scattering, and is propagate directly to reach and change
Energy device, and the signal being detected；So-called scattered wave refers to from from imaging region, by the one or many folding of scatterer
After penetrating, scattering, and propagates and reach energy converter, and the signal being detected；
The information vector of ultrasound transducer array is denoted as P (T, f)=[P1(T,f),…,Pn(T,f),…,PN(T, f)], directly
Arrived wave signal is denoted asScattered wave signal PS(T, f) is denoted as
The direct-path signal of n-th of transducer unitMeet:
Wherein, A0=Z1/2, k=2 π f/c, c are the medium velocity of sound；Z=cT indicates imaging plane to the vertical of transducer array
Distance also referred to as corresponds to the depth of time gated T；X indicates the x-axis coordinate value of arbitrary point on imaging plane, which indicates
For (X, Z)；
If there are L scattering paths between arbitrary point (X, Z) and n-th of transducer unit on imaging plane, then change for n-th
The scattered wave signal of energy device unitMeet:
Wherein, AlWith slFor the corresponding amplitude of the l articles scattering path and phase parameter；Since scatterer is random distribution, because
This AlWith slAlso corresponding randomness is embodied；
Step 32: information matrix K is constructed using information vector P:
Information matrix K is according to whether be divided into two parts containing random parameter: right side of the equal sign first item is denoted as matrix KC, etc.
Residual term on the right of number is denoted as matrix KR；
Matrix KCIt is without containing any random parameter, i.e., unrelated with the distributing position of scatterer, and have such as between matrix element
Lower correlation:
Wherein, m=1,2 ..., N, 1≤m ± q≤N；Representing matrix KCIn be located at back-diagonal and apart from diagonal line away from
From the phase difference between the element for q,Representing matrix KCIn be located at the position (m-q, m+q) at element,It indicates
Matrix KCIn be located at the position (m, m) at element；
Matrix KRMiddle all elements contain random parameter, show randomness, therefore do not have phase between matrix element
Guan Xing.
Scattered specifically, in the step 4, in information matrix K wave component to filter out operating process as follows:
Step 41: the method being different from document [1], the invention patent present invention is according to the following formula by the information matrix K inverse time
Needle rotates 45 °, and to respective element zero padding, to guarantee that matrix dimensionality is constant, obtains two homography A1={ a1u,v, A2=
{a2u,v}:
If N is odd number:
If N is even number:
Wherein, km,nIt indicates to be located at the element at the position (m, n) in information matrix K,
The dimension of matrix A 1
For N, the dimension of matrix A 2 is N -1；Matrix A 1 and the unification of matrix A 2 are indicated with matrix A, and use NAThe dimension of representing matrix A；
Step 42: construction filtering matrixMatrix A is filtered, the expression formula of filter vector C are as follows:
Wherein,For the conjugate transposition of C；U=1,2 ..., NA,
Step 43: matrix A being filtered, filtered matrix A is denoted as matrix AF:
AF=FA=F (AC+AR)=FAC+FAR
Wherein, matrix ACFor the part for being free of random parameter in matrix A, matrix ARFor in matrix A containing the portion of random parameter
Point；First item FAC=AC, matrix ACIt is remained unchanged before and after filtering, and matrix ARForcefully weakened after filtering is used, is dissipated
Ejected wave part is effectively filtered out；
Step 44: by matrix AF45 ° are rotated clockwise, filtered information matrix K is obtained, by filtered information matrix
K is denoted as matrix KF:
If (m-n) is odd number,
If (m-n) is even number,
Wherein,Representing matrix KFIn be located at the position (m, n) at element,Representing matrix A1F
In be located at the position ((m-n)/2+ (N+3)/4, (m+n)/2) at element,Representing matrix A2FMiddle position
Element at the position ((m-n-1)/2+ (N+3)/4, (m+n-1)/2), matrix A 1FFor filtered matrix A 1, matrix A 2FFor
Filtered matrix A 2；
Matrix KFSize be still N × N, it is consistent with information matrix K, filtering front and back without information dimension missing.
Specifically, it is as follows to carry out time reversal imaging process according to filtered information matrix in the step 5:
Step 51: to matrix KFMake singular value decomposition:
Wherein, UFFor N × N rank unitary matrice, its row vector is matrix KFLeft singular vector；ΛFFor positive semidefinite N × N rank
Diagonal matrix；For N × N rank unitary matrice,For VFConjugate transposition, VFRow vector be matrix KFRight singular vector；
ΛFElement on diagonal line is matrix KFSingular value；
Step 52: calculating the direct wave transmission matrix on imaging plane between arbitrary point (X, Z) and ultrasound transducer array
The element g being located at the position (m, n) in G, matrix Gm,nAre as follows:
Wherein, rm,nBetween n-th of the transducer unit and m-th of transducer unit in ultrasound transducer array away from
From；
Step 53: feature vector corresponding to maximum M singular value and the reversion of direct wave transmission matrix G application time are calculated
Son, i-th of singular valuePhotoacoustce signal intensity at depth Z are as follows:
Wherein,Indicate VFThe i-th row singular vector；The photoacoustic signal found out corresponding to maximum M singular value is strong
Degree superposition, obtains the photoacoustic image at depth Z It is IFM-th of element of vector, its table
Show preliminary at picture value with position corresponding to m-th of transducer unit on imaging plane.
Specifically, corrected parameter generates as follows with the process of utilization in the step 6:
Step 61: at depth Z, imaging correction value corresponding to m-th of transducer unitAre as follows:
Wherein,For gm,nConjugation；
Step 62: utilizing imaging correction valueTo tentatively at picture valueIt is modified, to obtain final at picture value
That is, final imaging results are by filtering imaging results IFWith correction matrix IRCorresponding element is divided by obtain.
The utility model has the advantages that the image weight of the Photoacoustic tomography provided by the invention based on matrix filtering and time reversal operator
Structure method has the advantages that the information matrix that the 1, present invention is received using ultrasound transducer array, can be effectively to sound wave
Signal is filtered, and is reduced interference of the scattering medium to imaging effect, is obtained the target imaging of high quality；2, the present invention only
Need to once excite and receive process, it is easy to operate efficiently, practicability with higher and wide applicability；3, the present invention is filtering
It not will cause the loss of information after wavefront, all effective informations received can be used；4, the present invention is according to uniformly
Acoustic transfer function in medium acquires corrected parameter, reduces pattern distortion, void caused by detecting due to limited perspective ultrasonic signal
The image deflects such as false contrast.
Detailed description of the invention
Fig. 1 is the principals of image reconstruction signal of the Optimization-type Photoacoustic tomography based on matrix filtering and time reversal operator
Figure；
Fig. 2 is matrix rotation, zero padding and filtering operation schematic diagram；
Fig. 3 is filtering front and back feature space matrix comparison diagram, 3 (a) the signal matrix K to receive in the case of no scattering layer
Real part image, 3 (b) the real part image for signal matrix K after scattering layer is added in a model, 3 (c) is using in document [2]
The filtered matrix K of filtering methodFReal part image, 3 (d) be using the filtered matrix of filtering method proposed by the invention
KFReal part image；
Fig. 4 is tradition filtering imaging method and Optimization-type filters imaging method imaging results comparison diagram, and 4 (a) be conventional wave
The imaging results of beam shaping imaging method, 4 (b) be the imaging results of time reversal in non-filtered situation, and 4 (c) be document [2]
The imaging results of middle filtering imaging method, 4 (d) filter the imaging results of imaging method for Optimization-type proposed by the invention.
Specific embodiment
The present invention will be further explained with reference to the accompanying drawing.
It is as shown in Figure 1 a kind of image reconstruction method of Photoacoustic tomography based on matrix filtering and time reversal operator
Schematic illustration includes the following steps for the small ball imaging after scattering layer.
Step 1: the small steel ball that 40 diameters are 0.8mm is arranged at random at random to be placed in as random scatter layer
Between ultrasonic transducer and target area.Due to the irregular distribution of scatterer, in the photoacoustic signal received, in addition to comprising coming
From in the direct wave of imageable target, also comprising via the scattering wave component after scattering layer Multiple Scattering.Scattering layer is dry to generate
Disturb signal.Four small steel balls that diameter is 0.8mm are placed in after scattering layer as imageable target.Quilt is irradiated using pulse laser
Imageable target after the small steel ball in imageable target absorbs laser energy, radiates ultrasonic wave to surrounding media due to optoacoustic effect
(system is as shown in Figure 1).
Step 2: the ultrasonic signal that sample emission goes out is received using the ultrasound transducer array that sample frequency is 20MHz, it should
Ultrasound transducer array has N=101 unit, and the spacing between adjacent cells is w=0.5mm, therefore the cross of n-th of unit
It is x to coordinaten=[n- (N+1)/2] w, N is odd number herein, and when N is even number, method of the invention is still set up.
Step 3: time gated and Short Time Fourier Transform being carried out to the ultrasonic signal received and is operated, it is obtained and corresponds to
The signal component of frequency f constructs the information matrix K of N × N.
Step 31: what interception ultrasonic transducer received is in time rangeInterior signal, and it is right
It makees Short Time Fourier Transform, obtains frequency domain presentation, is denoted as Pn(T, f), wherein subscript n indicates corresponding n-th of transducer unit.
Frequency-region signal configuration information vector P (T, the f)=[P obtained by N number of energy converter1(T,f),…,Pn(T,f),…,PN(T,f)]。
Information vector P (T, f) includes direct-path signal PD(T, f) and scattered wave signal PS(T, f) two parts (as shown in Figure 1).It is through
Wave signal corresponding element expression formula meets:
Wherein, A0=Z1/2, k=2 π f/c, c are the medium velocity of sound.Correspondingly, scattered wave signal corresponding element expression formula is full
Foot:
Wherein, AlWith slFor the corresponding amplitude of the l articles scattering path and phase parameter；Since scatterer is random distribution, because
This AlWith slAlso corresponding randomness is embodied.
Step 32: information matrix K is constructed using information vector P:
Matrix K is according to whether be divided into two parts containing random parameter.Right side of the equal sign first item is denoted as KC, this do not contain
Any random parameter, i.e., it is unrelated with the distributing position of scatterer, and there is following correlation between matrix element:
In contrast, random parameter is contained in the residual term of right side of the equal sign, system is denoted as KR。KRShow randomness, element
Between do not have KCThe correlation embodied.
Step 4: to information matrix K carry out 45 ° of rotations counterclockwise, zero padding, filtering, 45 ° rotate clockwise to restore etc. and operate,
Scattering wave component is filtered out, through wave component is retained.
Step 41: information matrix is rotated to 45 ° counterclockwise according to the following formula, obtains homography A1, A2:
M=u+v-(N+1)/2, n=v-u+ (N+1)/2 in above formula.The dimension of matrix A 1 is N, and the dimension of A2 is N -1.For
Narration is convenient, hereafter unifies to indicate A1 or A2 with A, and use NAThe dimension of representing matrix A.
Rotation, zero padding process as shown in Fig. 2, rotate 45 ° for original signal matrix counterclockwise.The coherence of original matrix is by opposing
Linea angulata direction is transferred to column direction.Then, postrotational matrix is divided by two parts, light green color circle according to the parity of columns
Shape part corresponds to odd column, and gray diamonds part corresponds to even column.Finally, two parts are separated, then respective zero padding is (light blue
Square portion), to obtain matrix A 1, A2.
Step 42: construction filtering matrixA matrix is filtered.The expression formula of filter vector C are as follows:
Wherein u=1,2 ..., NA,
Step 43: matrix A being filtered, filtering matrix A is obtainedF
AF=FA=F (AC+AR)=FAC+FAR
First item FA in above formulaC=AC, correlation item ACIt is remained unchanged before and after filtering, and random entry ARBy filtering square
Forcefully weakened after battle array effect, scattered wave part is effectively filtered out.
Step 44: by matrix AF45 ° are rotated clockwise, filtered information matrix K is obtainedF:
If (m-n) is odd number,
If (m-n) is even number,
The information matrix K obtained after filteringFSize be still N × N, it is consistent with original signal matrix K, filtering front and back do not have
The missing of information dimension.
When access time T=58 μ s (the small steel ball for corresponding to [- 5mm, 87mm]), centre frequency f=2.0MHz,
The contrast effect of matrix filtering front and back is as shown in Figure 3.In the case of Fig. 3 (a) is no scattering layer, the reality of the signal matrix K received
Portion's image.Matrix signal only contains direct wave ingredient at this time, and picture showing goes out clearly correlation.When we are added in a model
After scattering layer, receives the scattering wave component in signal and occupy leading position, shown in real part image such as Fig. 3 (b) of signal matrix K,
Image appearance is strong randomness at this time.Fig. 3 (c) is using the filtered matrix K of filtering method in document 1FReal part figure
Picture, the scatter-wave jamming of picture centre region obtained it is a degree of filter out, it can be observed that correlation pattern, however lack
The bulk informations of peripheral regions.Fig. 3 (d) is using the filtered matrix K of filtering method proposed by the inventionFReal part figure
Picture.It is observed that on the one hand the interference of scattering wave component has obtained more effectively filtering out, it can be observed that clearly correlation
Pattern, it is as a result highly close with 3 (a).On the other hand, matrix dimensionality remains unchanged, and the missing of information does not occur, after being
Imaging process provides better guarantee.
Step 5: the corresponding imaging for propagating section is obtained to filtered matrix application time reversal.
Step 51: to filtering matrix KFMake singular value decomposition:
Step 52: the direct wave transmission matrix G between the point and ultrasound transducer array unit on focal plane is calculated, yuan
Plain expression formula are as follows:
R in formulam,nIt is the distance between m-th point on n-th of unit of transducer array and focal plane, m, n=(1,
2,…,N)。
Step 53: feature vector corresponding to maximum preceding several singular values and the reversion of transmission matrix G application time are calculated
Son, i-th of singular valuePhotoacoustce signal intensity at depth Z are as follows:
The photoacoustce signal intensity found out corresponding to preceding M non-zero singular value is superimposed, the photoacoustic image of corresponding depth is obtained
Step 6: according to the acoustic transfer function in homogeneous medium, acquiring corrected parameter, reduce since limited perspective ultrasound is believed
The image deflects such as pattern distortion, false contrast caused by number detecting.
Step 61: at specified depth Z, imaging correction value corresponding to m-th point is obtained by following formula on focal plane:
Step 62: utilizing the correction valueTo tentatively at picture valueIt is modified, to obtain final at picture value
Final imaging results are by filtering imaging results IFWith correction matrix IRCorresponding element is divided by obtain:
Step 7: changing time gated region, repeat step 3~6, obtain the photoacoustic imaging of different depth, finally obtain whole
The photoacoustic image in a region.
Fig. 4 is the comparison diagram of a variety of imaging method imaging results, and Fig. 4 (a), 4 (b) be respectively traditional beam forming imaging method
And in non-filtered situation time reversal imaging results.As can be seen that the image quality of both methods is poor, background
There is a large amount of interference images and distortion in region, it is fairly obvious to scatter influence of the wave component to image quality.Fig. 4 (c) is in document 1
The contrast of the imaging results of filtering imaging method, absorber of light and background area significantly improves, and filtering is to image quality
Improvement emerge from.However this method only realizes imaging, two light absorptions of upper and lower two side areas to central area
Body fails to be imaged.Fig. 4 (d) is the imaging results that Optimization-type proposed by the invention filters imaging method.Due to effectively believing
The utilization rate of breath is improved, and the contrast of absorber of light and background area further increases, while obtaining complete area
Imaging results, all absorber of light have obtained effective imaging.
The above is only a preferred embodiment of the present invention, it should be pointed out that: for the ordinary skill people of the art
For member, various improvements and modifications may be made without departing from the principle of the present invention, these improvements and modifications are also answered
It is considered as protection scope of the present invention.
Claims (5)
1. a kind of image reconstructing method of the Photoacoustic tomography based on complete matrix filtering and time reversal operator, feature exist
In: include the following steps:
Step 1: using pulsed laser irradiation by Imaged samples, after the absorber of light in sample absorbs the energy of pulse laser, by
Ultrasonic wave is radiated to surrounding media in optoacoustic effect；
Step 2: receiving the ultrasonic signal that sample emission goes out using ultrasound transducer array, the ultrasound transducer array has
The transducer unit of N number of arrangement of straight line along the x-axis direction, the center spacing of two neighboring transducer unit are w, n-th of energy converter
The lateral coordinates of unit are as follows:
Step 3: the ultrasonic signal received is carried out including operation time gated, including Short Time Fourier Transform, construction
Size is the information matrix K of N × N；
Step 4: to information matrix K sequence carry out include 45 ° of rotations counterclockwise, zero padding, filtering, 45 ° rotate clockwise to restore and exist
Interior operation filters out scattering wave component, retains through wave component；
Step 5: singular value decomposition being carried out to filtered information matrix, to feature vector corresponding to maximum M singular value
Application time inverts operator, obtains the photoacoustic image of corresponding depth；
Step 6: according to the acoustic transfer function in homogeneous medium, acquiring corrected parameter, reduce since limited perspective ultrasonic signal is visited
Image deflects caused by surveying including pattern distortion, false contrast；
Step 7: changing time gated region, repeat step 3~6, obtain the photoacoustic imaging of different depth, finally obtain entire area
The photoacoustic image in domain.
2. the image reconstruction side of the Photoacoustic tomography described in claim 1 based on complete matrix filtering and time reversal operator
Method, it is characterised in that: in the step 3, the construction process of information matrix K is as follows:
Step 31: what interception ultrasonic transducer received is inSignal in time window, and it is made
Fourier transformation obtains the frequency domain presentation P of n-th of transducer unitn(T, f), frequency domain presentation Pn(T, f) believes comprising direct wave
Number PD(T, f) and scattered wave signal PSTwo parts (T, f)；Wherein, T indicates imaging plane to transducer array Acoustic Wave Propagation
It is shortest time, also referred to as time gated；Δ t indicates that the width of time window, f indicate frequency；
The information vector of ultrasound transducer array is denoted as P (T, f)=[P1(T,f),…,Pn(T,f),…,PN(T, f)], direct wave
Signal is denoted asScattered wave signal PS(T, f) is denoted as
The direct-path signal of n-th of transducer unitMeet:
Wherein, A0=Z1/2, k=2 π f/c, c are the medium velocity of sound；Z=cT indicate imaging plane to transducer array vertical range,
Also referred to as correspond to the depth of time gated T；X indicate imaging plane on arbitrary point x-axis coordinate value, the arbitrary point be expressed as (X,
Z)；
If there are L scattering path between arbitrary point (X, Z) and n-th of transducer unit on imaging plane, then n-th of energy converter
The scattered wave signal of unitMeet:
Wherein, AlWith slFor the corresponding amplitude of the l articles scattering path and phase parameter；
Step 32: information matrix K is constructed using information vector P:
Information matrix K is according to whether be divided into two parts containing random parameter: right side of the equal sign first item is denoted as matrix KC, the equal sign right side
The residual term on side is denoted as matrix KR；
Matrix KCIt is without containing any random parameter, i.e., unrelated with the distributing position of scatterer, and have between matrix element following related
Property:
Wherein, m=1,2 ..., N, 1≤m ± q≤N；Representing matrix KCIn be located at back-diagonal and apart from diagonal distance be q
Element between phase difference,Representing matrix KCIn be located at the position (m-q, m+q) at element,Representing matrix KC
In be located at the position (m, m) at element；
Matrix KRMiddle all elements contain random parameter, show randomness, therefore do not have correlation between matrix element.
3. the image reconstruction side of the Photoacoustic tomography described in claim 1 based on complete matrix filtering and time reversal operator
Method, it is characterised in that: scattered in the step 4, in information matrix K wave component to filter out operating process as follows:
Step 41: information matrix K being rotated 45 ° counterclockwise according to the following formula, and to respective element zero padding, to guarantee matrix dimensionality not
Become, obtains two homography A1={ a1u,v, A2={ a2u,v}:
If N is odd number:
If N is even number:
Wherein, km,nIt indicates to be located at the element at the position (m, n) in information matrix K,The dimension of matrix A 1 is N,
The dimension of matrix A 2 is N -1；Matrix A 1 and the unification of matrix A 2 are indicated with matrix A, and use NAThe dimension of representing matrix A；
Step 42: construction filtering matrixMatrix A is filtered, the expression formula of filter vector C are as follows:
Wherein,For the conjugate transposition of C；U=1,2 ..., NA,
Step 43: matrix A being filtered, filtered matrix A is denoted as matrix AF:
AF=FA=F (AC+AR)=FAC+FAR
Wherein, matrix ACFor the part for being free of random parameter in matrix A, matrix ARFor in matrix A containing the part of random parameter；The
One FAC=AC；
Step 44: by matrix AF45 ° are rotated clockwise, filtered information matrix K is obtained, filtered information matrix K is denoted as
Matrix KF:
If (m-n) is odd number,
If (m-n) is even number,
Wherein,Representing matrix KFIn be located at the position (m, n) at element,Representing matrix A1FMiddle position
Element at the position ((m-n)/2+ (N+3)/4, (m+n)/2),Representing matrix A2FIn be located at
Element at the position ((m-n-1)/2+ (N+3)/4, (m+n-1)/2), matrix A 1FFor filtered matrix A 1, matrix A 2FFor filter
Matrix A 2 after wave；
Matrix KFSize be still N × N, it is consistent with information matrix K, filtering front and back without information dimension missing.
4. the image reconstruction side of the Photoacoustic tomography described in claim 1 based on complete matrix filtering and time reversal operator
Method, it is characterised in that: in the step 5, it is as follows that time reversal imaging process is carried out according to filtered information matrix:
Step 51: to matrix KFMake singular value decomposition:
Wherein, UFFor N × N rank unitary matrice, its row vector is matrix KFLeft singular vector；ΛFIt is diagonal for positive semidefinite N × N rank
Matrix；For N × N rank unitary matrice,For VFConjugate transposition, VFRow vector be matrix KFRight singular vector；ΛFIt is right
Element on linea angulata is matrix KFSingular value；
Step 52: calculating the direct wave transmission matrix G on imaging plane between arbitrary point (X, Z) and ultrasound transducer array, square
The element g being located at the position (m, n) in battle array Gm,nAre as follows:
Wherein, rm,nFor n-th of the transducer unit and the distance between m-th of transducer unit in ultrasound transducer array；
Step 53: feature vector corresponding to maximum M singular value and direct wave transmission matrix G application time invert operator,
I-th of singular valuePhotoacoustce signal intensity at depth Z are as follows:
Wherein, Vi FIndicate VFThe i-th row singular vector；The photoacoustce signal intensity found out corresponding to maximum M singular value is folded
Add, obtains the photoacoustic image at depth Z It is IFM-th of element of vector, it is expressed as
As preliminary at picture value with position corresponding to m-th of transducer unit in plane.
5. the image reconstruction side of the Photoacoustic tomography described in claim 1 based on complete matrix filtering and time reversal operator
Method, it is characterised in that: in the step 6, corrected parameter generates as follows with the process of utilization:
Step 61: at depth Z, imaging correction value corresponding to m-th of transducer unitAre as follows:
Wherein,For gm,nConjugation；
Step 62: utilizing imaging correction valueTo tentatively at picture valueIt is modified, to obtain final at picture value
That is, final imaging results are by filtering imaging results IFWith correction matrix IRCorresponding element is divided by obtain.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810706816.4A CN109091109B (en) | 2018-07-02 | 2018-07-02 | Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810706816.4A CN109091109B (en) | 2018-07-02 | 2018-07-02 | Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109091109A true CN109091109A (en) | 2018-12-28 |
| CN109091109B CN109091109B (en) | 2021-04-20 |
Family
ID=64845298
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810706816.4A Active CN109091109B (en) | 2018-07-02 | 2018-07-02 | Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109091109B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109674490A (en) * | 2019-01-17 | 2019-04-26 | 南京大学深圳研究院 | A kind of low reflection artifacts photoacoustic microscope imaging method of ultrasonic guidance |
| CN109864707A (en) * | 2019-01-17 | 2019-06-11 | 南京科技职业学院 | A method of improving Photoacoustic tomography resolution ratio in limited perspective |
| CN110881947A (en) * | 2019-12-06 | 2020-03-17 | 北京信息科技大学 | Optical coherence tomography imaging method |
| CN111435528A (en) * | 2019-01-15 | 2020-07-21 | 中国科学院金属研究所 | Laser ultrasonic visual image quality improvement processing method |
Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040059265A1 (en) * | 2002-09-12 | 2004-03-25 | The Regents Of The University Of California | Dynamic acoustic focusing utilizing time reversal |
| US20040128081A1 (en) * | 2002-12-18 | 2004-07-01 | Herschel Rabitz | Quantum dynamic discriminator for molecular agents |
| CN101173869A (en) * | 2006-05-05 | 2008-05-07 | 尤洛考普特公司 | Method and apparatus for diagnosing a mechanism |
| CN101502131A (en) * | 2006-08-10 | 2009-08-05 | 皇家飞利浦电子股份有限公司 | A device for and a method of processing an audio signal |
| CN102068277A (en) * | 2010-12-14 | 2011-05-25 | 哈尔滨工业大学 | Method and device for observing photoacoustic imaging in single-array element and multi-angle mode based on compressive sensing |
| CN102132149A (en) * | 2008-06-27 | 2011-07-20 | 沃尔弗拉姆·R·雅里施 | High efficiency computed tomography |
| CN102867294A (en) * | 2012-05-28 | 2013-01-09 | 天津大学 | Fourier-wavelet regularization-based coaxial phase contrast image restoration method |
| US20150249493A1 (en) * | 2012-09-24 | 2015-09-03 | Telefonaktiebolaget L M Ericsson (Publ) | Prefiltering in MIMO Receiver |
| CN105181805A (en) * | 2015-09-30 | 2015-12-23 | 中国计量学院 | Multi-filtering ultrasonic imaging method based on time-reversal MUSIC |
| CN105785566A (en) * | 2016-03-31 | 2016-07-20 | 南京大学 | Method utilizing spatial light modulator to improve photoacoustic imaging limited view angle |
| CN107710014A (en) * | 2015-05-12 | 2018-02-16 | 法国电力公司 | The method and apparatus detected is propagated using ripple |
| CN107861517A (en) * | 2017-11-01 | 2018-03-30 | 北京航空航天大学 | The online trajectory planning method of guidance of great-jump-forward reentry vehicle based on linear pseudo- spectrum |
-
2018
- 2018-07-02 CN CN201810706816.4A patent/CN109091109B/en active Active
Patent Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040059265A1 (en) * | 2002-09-12 | 2004-03-25 | The Regents Of The University Of California | Dynamic acoustic focusing utilizing time reversal |
| US20040128081A1 (en) * | 2002-12-18 | 2004-07-01 | Herschel Rabitz | Quantum dynamic discriminator for molecular agents |
| CN101173869A (en) * | 2006-05-05 | 2008-05-07 | 尤洛考普特公司 | Method and apparatus for diagnosing a mechanism |
| CN101502131A (en) * | 2006-08-10 | 2009-08-05 | 皇家飞利浦电子股份有限公司 | A device for and a method of processing an audio signal |
| CN102132149A (en) * | 2008-06-27 | 2011-07-20 | 沃尔弗拉姆·R·雅里施 | High efficiency computed tomography |
| CN102068277A (en) * | 2010-12-14 | 2011-05-25 | 哈尔滨工业大学 | Method and device for observing photoacoustic imaging in single-array element and multi-angle mode based on compressive sensing |
| CN102867294A (en) * | 2012-05-28 | 2013-01-09 | 天津大学 | Fourier-wavelet regularization-based coaxial phase contrast image restoration method |
| US20150249493A1 (en) * | 2012-09-24 | 2015-09-03 | Telefonaktiebolaget L M Ericsson (Publ) | Prefiltering in MIMO Receiver |
| CN107710014A (en) * | 2015-05-12 | 2018-02-16 | 法国电力公司 | The method and apparatus detected is propagated using ripple |
| CN105181805A (en) * | 2015-09-30 | 2015-12-23 | 中国计量学院 | Multi-filtering ultrasonic imaging method based on time-reversal MUSIC |
| CN105785566A (en) * | 2016-03-31 | 2016-07-20 | 南京大学 | Method utilizing spatial light modulator to improve photoacoustic imaging limited view angle |
| CN107861517A (en) * | 2017-11-01 | 2018-03-30 | 北京航空航天大学 | The online trajectory planning method of guidance of great-jump-forward reentry vehicle based on linear pseudo- spectrum |
Non-Patent Citations (3)
| Title |
|---|
| CHAO HUANG等: "Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography With Acoustically Inhomogeneous Media", 《IEEE TRANSACTIONS ON MEDICAL IMAGING》 * |
| WEI RUI等: "Photoacoustic imaging in scattering media by combining a correlation matrix filter with a time reversal operator", 《OPTICS EXPRESS》 * |
| 姜自波 等: "快速在体光声计算层析图像重建方法", 《计算机应用》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111435528A (en) * | 2019-01-15 | 2020-07-21 | 中国科学院金属研究所 | Laser ultrasonic visual image quality improvement processing method |
| CN109674490A (en) * | 2019-01-17 | 2019-04-26 | 南京大学深圳研究院 | A kind of low reflection artifacts photoacoustic microscope imaging method of ultrasonic guidance |
| CN109864707A (en) * | 2019-01-17 | 2019-06-11 | 南京科技职业学院 | A method of improving Photoacoustic tomography resolution ratio in limited perspective |
| CN109864707B (en) * | 2019-01-17 | 2021-09-07 | 南京科技职业学院 | Method for improving photoacoustic tomography resolution ratio under limited viewing angle condition |
| CN109674490B (en) * | 2019-01-17 | 2021-09-10 | 南京大学深圳研究院 | Ultrasonic-guided photoacoustic microscope imaging method with low reflection artifact |
| CN110881947A (en) * | 2019-12-06 | 2020-03-17 | 北京信息科技大学 | Optical coherence tomography imaging method |
| CN110881947B (en) * | 2019-12-06 | 2022-07-05 | 北京信息科技大学 | Optical coherence tomography imaging method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109091109B (en) | 2021-04-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109091109A (en) | The image reconstructing method of Optimization-type Photoacoustic tomography based on complete matrix filtering and time reversal operator | |
| OˈReilly et al. | A super‐resolution ultrasound method for brain vascular mapping | |
| Gateau et al. | Three‐dimensional optoacoustic tomography using a conventional ultrasound linear detector array: Whole‐body tomographic system for small animals | |
| Wang | Prospects of photoacoustic tomography | |
| US10653321B2 (en) | Photoacoustic computed tomography with a biplanar acoustic reflector | |
| US20040059265A1 (en) | Dynamic acoustic focusing utilizing time reversal | |
| Song et al. | Coded excitation plane wave imaging for shear wave motion detection | |
| Abadi et al. | Frequency-sum beamforming for passive cavitation imaging | |
| Ding et al. | Ultrasound line-by-line scanning method of spatial–temporal active cavitation mapping for high-intensity focused ultrasound | |
| Hemmsen et al. | Tissue harmonic synthetic aperture ultrasound imaging | |
| Lu et al. | Two-step aberration correction: application to transcranial histotripsy | |
| RU2378989C2 (en) | Method of diagnostics by means of ultrasonic, sonic and electromagnetic waves | |
| Francis et al. | Multiview spatial compounding using lens-based photoacoustic imaging system | |
| Gao et al. | Achieving depth-independent lateral resolution in AR-PAM using the synthetic-aperture focusing technique | |
| Warbal et al. | In silico evaluation of the effect of sensor directivity on photoacoustic tomography imaging | |
| Mao et al. | Improving photoacoustic imaging in low signal-to-noise ratio by using spatial and polarity coherence | |
| Lu et al. | Passive cavitation mapping using dual apodization with cross-correlation in ultrasound therapy monitoring | |
| Nabavizadeh et al. | In vivo Young's modulus mapping of pancreatic ductal adenocarcinoma during HIFU ablation using harmonic motion elastography (HME) | |
| Zheng et al. | Volumetric tri-modal imaging with combined photoacoustic, ultrasound, and shear wave elastography | |
| US12016657B2 (en) | Devices and methods for photoacoustic tomography | |
| Dean-Ben et al. | High resolution transcranial imaging based on the optoacoustic memory effect | |
| Kurnikov et al. | Fisheye piezo polymer detector for scanning optoacoustic angiography of experimental neoplasms | |
| Lu et al. | Dual apodization with cross‐correlation combined with robust Capon beamformer applied to ultrasound passive cavitation mapping | |
| Ahmed et al. | Parallel receive beamforming improves the performance of focused transmit-based single-track location shear wave elastography | |
| Paul et al. | Improvement of LED‐based photoacoustic imaging using lag‐coherence factor (LCF) beamforming |
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 |