阅读说明：本技术 一种生物腔体光声内窥成像方法及系统 (Photoacoustic endoscopic imaging method and system for biological cavity ) 是由 孙正 段爽 于 2019-08-29 设计创作，主要内容包括：本发明公开一种生物腔体光声内窥成像方法及系统,获取生物腔体的位置r处的光吸收能量初始值A_0(r)和声速初始值c_(s,0)(r),分别计算位置r处的第k次的光吸收能量A_k(r)、第k+1次的光吸收能量A_(k+1)(r)、第k次的声速c_(s,k)(r)和第k+1次的声速c_(s,k+1)(r)；通过判断光吸收能量绝对差ε_(A,k)和声速绝对差是否满足收敛条件,决定是继续迭代还是直接构建光吸收能量分布图和声速分布图。本发明考虑生物组织的复杂性,综合考虑速度的差异性,结合光吸收能量和不同声速构建光声图像,减少由待测组织内的声速分布不均匀导致的误差,得到高质量的腔体横截面光吸收能量分布图。(The invention discloses a photoacoustic endoscopic imaging method and a photoacoustic endoscopic imaging system for a biological cavity, which are used for acquiring an initial value A of light absorption energy at a position r of the biological cavity 0 (r) and initial value of sound velocity c s,0 (r) calculating light absorption energy A of kth time at position r k (r) and (k + 1) th light absorption energy A k+1 (r), speed of sound c of kth s,k (r) and sound velocity c of (k + 1) th order s,k+1 (r); by judging the absolute difference epsilon of light absorption energy A,k And absolute difference of sound velocity And whether a convergence condition is met or not is determined, and whether a light absorption energy distribution graph and an acoustic velocity distribution graph are continuously iterated or directly constructed. The invention considers the complexity of biological tissues, comprehensively considers the difference of the speeds, combines light absorption energy and different sound velocities to construct a photoacoustic image, and reduces the problem that the sound velocity in the tissue to be detected is not distributedAnd errors caused by uniformity are avoided, and a high-quality cavity cross section light absorption energy distribution diagram is obtained.)

1. A photoacoustic endoscopic imaging method for a biological cavity, the method comprising:

obtaining the initial value A of the light absorption energy at the position r of the biological cavity₀(r) and initial value of sound velocity c_s,0(r)；

Respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r); k is iteration times, and k is more than or equal to 1;

respectively calculating the sound velocity c of the kth time at the position r of the biological cavity_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r)；

Calculating the light absorption energy A of the k time_k(r) and the light absorption energy A of the k +1 th order_k+1Absolute difference of light absorption energy ε of (r)_A,kAnd the speed of sound c of the kth time_s,k(r) speed of sound c at the k +1 th time_s,k+1Absolute difference in sound velocity of (r)

Judging the absolute difference epsilon of the light absorption energy_A,kAnd the absolute difference of sound velocityWhether the convergence condition is met or not is judged, and a convergence judgment result is obtained;

if the convergence judgment result is yes, outputting the (k + 1) th light absorption energy A at the position r_k+1(r) and speed of sound c_s,k+1(r) respectively constructing a light absorption energy distribution map and an acoustic velocity distribution map;

if the convergence judgment result is negative, updating the iteration times, and returning to the step of respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r)。

2. The photoacoustic endoscopic imaging method according to claim 1, wherein the formula A is given as the formula_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r)) calculating the light absorption energy A of the k-th time_k(r) and the light absorption energy A of the k +1 th order_k+1(r)；

Wherein the content of the first and second substances,

is a first photoacoustic relationship function;

arg min[·]is that g (A)_k(r))+λΦ_TV(A_k(r)) A at the minimum_k(r)；

W(A_k-1(r)) is an approximate blackplug matrix of the first photoacoustic relationship function;

W^-1(A_k-1(r)) is W (A)_k-1(r)) an inverse matrix;

a is a light absorption energy distribution matrix; a is more than or equal to 0, which means the light absorption energy which is more than or equal to 0 in the light absorption energy distribution matrix;

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²is a second photoacoustic relationship function;

p_m(r) is a measurement of the photoacoustic signal at location r;

H(c_s,k-1(r)) is a first operator related to the speed of sound;

| | · | is a 2-norm;

Φ_TV(A_k(r)) is a TV regularization term,

λ is the TV regularization parameter;

L_kis about A_k(r) a sparse matrix of variance characteristics;

η > 0 is a constant;

G₁'(r,A_k-1(r),c_s,k-1(r)) is the gradient of the photoacoustic relationship function.

3. Biological chamber according to claim 2The body photoacoustic endoscopic imaging method is characterized by being based on a formulaRespectively calculating the sound velocity c of the k-th time_s,k(r) and the sound speed c of the k +1 st time_s,k+1(r)；

Wherein the content of the first and second substances,

is f (q)_k) A gradient of (a);

f(q_k)＝||p_m(r)-H(q_k)·A_k-1(r)||²；

γ_k-1is the weight parameter, γ, after the k-1 iteration_k-2Is the weight parameter after the k-2 iteration;

weight parameter gamma after the kth iteration_kIs composed of

H(q_k) Is a second operator related to the speed of sound;

d is f (c)_s,k-1(r)) the prestz constant of the derivative;

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²。

4. the photoacoustic endoscopic imaging method according to claim 1, wherein the convergence condition is ε_A,k＜ε_AAnd is

Wherein epsilon_AFor the convergence tolerance of the light absorption energy,is the convergence tolerance of the speed of sound.

5. The photoacoustic endoscopic imaging method according to claim 1, wherein the constructing the light absorption energy distribution image and the sound velocity distribution image specifically comprises:

the light absorption energy distribution matrix A under the polar coordinates is normalized and subjected to gray scale processing by adopting the following formula:

b (i, j) is a value obtained by normalizing and graying a (i, j);

i is the abscissa of the biological cavity position r under the corresponding polar coordinate system, and j is the ordinate of the biological cavity position r under the corresponding polar coordinate system;

a (i, j) is the element value of the ith row and the jth column in the light absorption energy distribution matrix A;

min (a) is the minimum value of the elements in the light absorption energy distribution matrix a, max (a) is the maximum value of the elements in the light absorption energy distribution matrix a;

adopting the following formula to align the sound velocity distribution matrix c under polar coordinates_sCarrying out normalization and gray level processing;

e (i, j) is for c_s(i, j) normalized and grayed values;

c_s(i, j) is the sound velocity distribution matrix c_sRow (i) and column (j)The value of the element(s);

min(c_s) Is the sound velocity distribution matrix c_sMinimum value of (1); max (c)_s) Is the sound velocity distribution matrix c_sMaximum value of (1);

converting the B (i, j) into B '(x, y) in plane rectangular coordinates according to the formula B' (x, y) ═ B (jcosi, jsini);

converting the E (i, j) into E '(x, y) in plane rectangular coordinates according to the formula E' (x, y) ═ E (jcosi, jsini);

wherein i belongs to [0,360], j belongs to [0, d ];

d is the maximum value of the polar diameter in the theta-l polar coordinate system;

b' (x, y) is a gray value of a point (x, y) in a rectangular coordinate system in the light absorption energy distribution map; e' (x, y) is the grayscale value of the point (x, y) in the rectangular coordinate system in the sound velocity map.

6. A biological cavity photoacoustic endoscopic imaging system, the system comprising:

an initial value acquisition module for acquiring an initial value A of light absorption energy at a position r of the biological cavity₀(r) and initial value of sound velocity c_s,0(r)；

A light absorption energy calculation module for respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r); k is iteration times, and k is more than or equal to 1;

a sound velocity calculation module for respectively calculating the kth sound velocity c at the position r of the biological cavity_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r)；

An absolute difference calculation module for calculating the light absorption energy A of the kth time_k(r) and the light absorption energy A of the k +1 th order_k+1Absolute difference of light absorption energy ε of (r)_A,kAnd the speed of sound c of the kth time_s,k(r) speed of sound c at the k +1 th time_s,k+1Absolute difference in sound velocity of (r)

A convergence judgment result module for judging the absolute difference epsilon of the light absorption energy_A,kAnd the absolute difference of sound velocityWhether the convergence condition is met or not is judged, and a convergence judgment result is obtained;

an image construction module connected with the convergence judgment result module and used for absorbing energy A according to the (k + 1) th light at the position r when the convergence judgment result module is yes_k+1(r) and speed of sound c_s,k+1(r) respectively constructing a light absorption energy distribution map and an acoustic velocity distribution map;

and the iteration updating module is respectively connected with the convergence judgment result module and the light absorption energy calculating module, and is used for updating the iteration times and returning to the light absorption energy calculating module when the convergence judgment result module is negative.

7. The bio-cavity photoacoustic endoscopic imaging system according to claim 6, wherein formula A is given as the basis_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r)) calculating the light absorption energy A of the k-th time_k(r) and the light absorption energy A of the k +1 th order_k+1(r)；

Wherein the content of the first and second substances,

is a first photoacoustic relationship function;

arg min[·]is that g (A)_k(r))+λΦ_TV(A_k(r)) A at the minimum_k(r)；

W(A_k-1(r)) is an approximate blackplug matrix of the first photoacoustic relationship function;

W^-1(A_k-1(r)) is W (A)_k-1(r)) an inverse matrix;

a is a light absorption energy distribution matrix; a is more than or equal to 0, which means the light absorption energy which is more than or equal to 0 in the light absorption energy distribution matrix;

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²is a second photoacoustic relationship function;

p_m(r) is a measurement of the photoacoustic signal at location r;

H(c_s,k-1(r)) is a first operator related to the speed of sound;

| | · | is a 2-norm;

Φ_TV(A_k(r)) is a TV regularization term,

λ is the TV regularization parameter;

L_kis about A_k(r) a sparse matrix of variance characteristics;

η > 0 is a constant;

G₁'(r,A_k-1(r),c_s,k-1(r)) is the gradient of the photoacoustic relationship function.

8. The bio-cavity photoacoustic endoscopic imaging system according to claim 7, wherein the formula is givenRespectively calculating the sound velocity c of the k-th time_s,k(r) and the speed of sound c after the (k + 1) th iteration_s,k+1(r)；

Wherein the content of the first and second substances,

is f (q)_k) A gradient of (a);

f(q_k)＝||p_m(r)-H(q_k)·A_k-1(r)||²；

γ_k-1is the weight parameter, γ, after the k-1 iteration_k-2Is the weight parameter after the k-2 iteration;

weight parameter gamma after the kth iteration_kIs composed of

H(q_k) Is a second operator related to the speed of sound;

d is f (c)_s,k-1(r)) the prestz constant of the derivative;

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²。

9. the photoacoustic endoscopic imaging system of claim 6 wherein the convergence condition is ε_A,k＜ε_AAnd is

Wherein epsilon_AFor the convergence tolerance of the light absorption energy,is the convergence tolerance of the speed of sound.

10. The bio-cavity photoacoustic endoscopic imaging system according to claim 6, wherein the image construction module specifically comprises:

the light absorption energy distribution matrix processing unit is used for carrying out normalization and gray scale processing on the light absorption energy distribution matrix A under the polar coordinates by adopting the following formula:

b (i, j) is a value obtained by normalizing and graying a (i, j);

i is the abscissa of the biological cavity position r under the corresponding polar coordinate system, and j is the ordinate of the biological cavity position r under the corresponding polar coordinate system;

a (i, j) is the element value of the ith row and the jth column in the light absorption energy distribution matrix A;

min (a) is the minimum value of the elements in the light absorption energy distribution matrix a, max (a) is the maximum value of the elements in the light absorption energy distribution matrix a;

a sound velocity distribution matrix processing unit for aligning the sound velocity distribution matrix c under polar coordinates by using the following formula_sCarrying out normalization and gray level processing;

e (i, j) is for c_s(i, j) normalized and grayed values;

c_s(i, j) is the sound velocity distribution matrix c_sThe value of the element in the ith row and the jth column;

min(c_s) Is the sound velocity distribution matrix c_sMinimum value of (1); max (c)_s) Is the sound velocity distribution matrix c_sMaximum value of (1);

a light absorption energy distribution map gray value conversion unit, which is used for converting the B (i, j) into B '(x, y) in a plane rectangular coordinate according to a formula B' (x, y) ═ B (jcosi, jsini);

a sound velocity profile gray value conversion unit for converting the E (i, j) into E '(x, y) in a plane rectangular coordinate according to the formula E' (x, y) ═ E (jcosi, jsini);

wherein i belongs to [0,360], j belongs to [0, d ];

d is the maximum value of the polar diameter in the theta-l polar coordinate system;

b' (x, y) is a gray value of a point (x, y) in a rectangular coordinate system in the light absorption energy distribution map; e' (x, y) is the grayscale value of the point (x, y) in the rectangular coordinate system in the sound velocity map.

Technical Field

The invention relates to the technical field of medical imaging, in particular to a photoacoustic endoscopic imaging method and a photoacoustic endoscopic imaging system for a biological cavity.

Background

A biological Photoacoustic endoscopic (PAE) imaging method is a novel medical functional imaging method, and is characterized in that a specially-made imaging catheter is directly inserted into a target cavity, a catheter probe emits short-pulse laser light to irradiate surrounding tissues, the tissues are heated and expanded after absorbing light energy, ultrasonic waves are excited and transmitted to the surfaces of the tissues, and Photoacoustic signals are obtained. An ultrasound probe at the tip of the imaging catheter receives photoacoustic signals from various directions during the intraluminal scan of the surrounding tissue. After the photoacoustic pressure time sequences collected at different measurement positions are sent to a computer, a light absorption energy distribution image of the cross section of the cavity can be inverted and reconstructed by adopting a proper image reconstruction algorithm, and the morphological structure and the functional components of the tissue are reflected.

Given the complexity of biological tissue, there is often a large difference in the speed of ultrasound as it propagates through tissue having different compositions. In many cases, the sound velocity distribution in the tissue cannot be predicted before PAE imaging is performed, so that the assumption of a constant sound velocity in the process of reconstructing the image may cause serious problems such as acoustic distortion, artifacts, blurring, and object dislocation in the reconstructed image.

Disclosure of Invention

The invention aims to provide a photoacoustic endoscopic imaging method and a photoacoustic endoscopic imaging system for a biological cavity, which are used for reducing errors caused by uneven sound velocity distribution in tissues to be measured and obtaining a high-quality light absorption energy distribution image of the cross section of the cavity.

In order to achieve the purpose, the invention provides the following scheme:

a method of photoacoustic endoscopic imaging of a biological cavity, the method comprising:

obtaining the initial value A of the light absorption energy at the position r of the biological cavity₀(r) and initial value of sound velocity c_s,0(r)；

Respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r); k is iteration times, and k is more than or equal to 1;

respectively calculating the sound velocity c of the kth time at the position r of the biological cavity_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r)；

Calculating the light absorption energy A of the k time_k(r) and the light absorption energy A of the k +1 th order_k+1Absolute difference of light absorption energy ε of (r)_A,kAnd the speed of sound c of the kth time_s,k(r) speed of sound c at the k +1 th time_s,k+1Absolute difference in sound velocity of (r)

Judging the absolute difference epsilon of the light absorption energy_A,kAnd the absolute difference of sound velocityWhether the convergence condition is met or not is judged, and a convergence judgment result is obtained;

if the convergence judgment result is yes, outputting the (k + 1) th light absorption energy A at the position r_k+1(r) and speed of sound c_s,k+1(r) respectively constructing a light absorption energy distribution map and an acoustic velocity distribution map;

if the convergence judgment result is negative, updating the iteration times, and returning to the step of respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r)。

Optionally, according to formula A_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r)) calculating the light absorption energy A of the k-th time_k(r) and the light absorption energy A of the k +1 th order_k+1(r)；

Wherein the content of the first and second substances,

argmin[·]is that g (A)_k(r))+λΦ_TV(A_k(r)) A at the minimum_k(r)；

W(A_k-1(r)) is an approximate blackplug matrix of the first photoacoustic relationship function;

W^-1(A_k-1(r)) is W (A)_k-1(r)) an inverse matrix;

a is a light absorption energy distribution matrix; a is more than or equal to 0, which means the light absorption energy which is more than or equal to 0 in the light absorption energy distribution matrix;

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²is a second photoacoustic relationship function;

p_m(r) is a measurement of the photoacoustic signal at location r;

H(c_s,k-1(r)) is a first operator related to the speed of sound;

| | · | is a 2-norm;

Φ_TV(A_k(r)) is a TV regularization term,

λ is the TV regularization parameter;

L_kis about A_k(r) a sparse matrix of variance characteristics;

η > 0 is a constant;

G₁'(r,A_k-1(r),c_s,k-1(r)) is the gradient of the photoacoustic relationship function.

Optionally, according to a formula

Respectively calculating the sound velocity c of the k-th time_s,k(r) and the sound speed c of the k +1 st time_s,k+1(r)；

Wherein the content of the first and second substances,

▽f(q_k) Is f (q)_k) A gradient of (a);

f(q_k)＝||p_m(r)-H(q_k)·A_k-1(r)||²；

γ_k-1is the weight parameter, γ, after the k-1 iteration_k-2Is the kth-Weight parameters after 2 iterations;

weight parameter gamma after the kth iteration_kIs composed of

H(q_k) Is a second operator related to the speed of sound;

d is f (c)_s,k-1(r)) the prestz constant of the derivative;

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²。

optionally, the convergence condition is epsilon_A,k＜ε_AAnd is

Wherein epsilon_AFor the convergence tolerance of the light absorption energy,is the convergence tolerance of the speed of sound.

Optionally, the constructing the light absorption energy distribution image and the sound velocity distribution image specifically includes:

the light absorption energy distribution matrix A under the polar coordinates is normalized and subjected to gray scale processing by adopting the following formula:

b (i, j) is a value obtained by normalizing and graying a (i, j);

i is the abscissa of the biological cavity position r under the corresponding polar coordinate system, and j is the ordinate of the biological cavity position r under the corresponding polar coordinate system;

a (i, j) is the element value of the ith row and the jth column in the light absorption energy distribution matrix A;

min (a) is the minimum value of the elements in the light absorption energy distribution matrix a, max (a) is the maximum value of the elements in the light absorption energy distribution matrix a;

adopting the following formula to align the sound velocity distribution matrix c under polar coordinates_sCarrying out normalization and gray level processing;

e (i, j) is for c_s(i, j) normalized and grayed values;

c_s(i, j) is the sound velocity distribution matrix c_sThe value of the element in the ith row and the jth column;

min(c_s) Is the sound velocity distribution matrix c_sMinimum value of (1); max (c)_s) Is the sound velocity distribution matrix c_sMaximum value of (1);

converting the B (i, j) into B '(x, y) in plane rectangular coordinates according to the formula B' (x, y) ═ B (jcosi, jsini);

converting the E (i, j) into E '(x, y) in plane rectangular coordinates according to the formula E' (x, y) ═ E (jcosi, jsini);

wherein i belongs to [0,360], j belongs to [0, d ];

d is the maximum value of the polar diameter in the theta-l polar coordinate system;

b' (x, y) is a gray value of a point (x, y) in a rectangular coordinate system in the light absorption energy distribution map; e' (x, y) is the grayscale value of the point (x, y) in the rectangular coordinate system in the sound velocity map.

A bio-cavity photoacoustic endoscopic imaging system, the system comprising:

an initial value acquisition module for acquiring an initial value A of light absorption energy at a position r of the biological cavity₀(r) and initial value of sound velocity c_s,0(r)；

A light absorption energy calculation module for respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r); k is iteration times, and k is more than or equal to 1;

a sound velocity calculation module for respectively calculating the kth sound velocity c at the position r of the biological cavity_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r)；

An absolute difference calculation module for calculating the light absorption energy A of the kth time_k(r) and the light absorption energy A of the k +1 th order_k+1Absolute difference of light absorption energy ε of (r)_A,kAnd the speed of sound c of the kth time_s,k(r) speed of sound c at the k +1 th time_s,k+1Absolute difference in sound velocity of (r)

A convergence judgment result module for judging the absolute difference epsilon of the light absorption energy_A,kAnd the absolute difference of sound velocityWhether the convergence condition is met or not is judged, and a convergence judgment result is obtained;

an image construction module connected with the convergence judgment result module and used for absorbing energy A according to the (k + 1) th light at the position r when the convergence judgment result module is yes_k+1(r) and speed of sound c_s,k+1(r) respectively constructing a light absorption energy distribution map and an acoustic velocity distribution map;

and the iteration updating module is respectively connected with the convergence judgment result module and the light absorption energy calculating module, and is used for updating the iteration times and returning to the light absorption energy calculating module when the convergence judgment result module is negative.

Optionally, according to formula A_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r)) calculating the light absorption energy A of the k-th time_k(r) and the light absorption energy A of the k +1 th order_k+1(r)；

Wherein the content of the first and second substances,

argmin[·]is that g (A)_k(r))+λΦ_TV(A_k(r)) A at the minimum_k(r)；

W(A_k-1(r)) is an approximate blackplug matrix of the first photoacoustic relationship function;

W^-1(A_k-1(r)) is W (A)_k-1(r)) an inverse matrix;

a is a light absorption energy distribution matrix; a is more than or equal to 0, which means the light absorption energy which is more than or equal to 0 in the light absorption energy distribution matrix;

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²is a second photoacoustic relationship function;

p_m(r) is a measurement of the photoacoustic signal at location r;

H(c_s,k-1(r)) is a first operator related to the speed of sound;

| | · | is a 2-norm;

Φ_TV(A_k(r)) is a TV regularization term,

λ is the TV regularization parameter;

L_kis about A_k(r) a sparse matrix of variance characteristics;

η > 0 is a constant;

G₁'(r,A_k-1(r),c_s,k-1(r)) is the gradient of the photoacoustic relationship function.

Optionally, according to a formula

Respectively calculating the sound velocity c of the k-th time_s,k(r) and the speed of sound c after the (k + 1) th iteration_s,k+1(r)；

Wherein the content of the first and second substances,

▽f(q_k) Is f (q)_k) A gradient of (a);

f(q_k)＝||p_m(r)-H(q_k)·A_k-1(r)||²；

γ_k-1is the weight parameter, γ, after the k-1 iteration_k-2Is the weight parameter after the k-2 iteration;

weight parameter gamma after the kth iteration_kIs composed of

H(q_k) Is a second operator related to the speed of sound;

d is f (c)_s,k-1(r)) the prestz constant of the derivative;

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²。

optionally, the convergence condition is epsilon_A,k＜ε_AAnd is

Wherein epsilon_AFor the convergence tolerance of the light absorption energy,is the convergence tolerance of the speed of sound.

Optionally, the image construction module specifically includes:

the light absorption energy distribution matrix processing unit is used for carrying out normalization and gray scale processing on the light absorption energy distribution matrix A under the polar coordinates by adopting the following formula:

b (i, j) is a value obtained by normalizing and graying a (i, j);

i is the abscissa of the biological cavity position r under the corresponding polar coordinate system, and j is the ordinate of the biological cavity position r under the corresponding polar coordinate system;

a (i, j) is the element value of the ith row and the jth column in the light absorption energy distribution matrix A;

min (a) is the minimum value of the elements in the light absorption energy distribution matrix a, max (a) is the maximum value of the elements in the light absorption energy distribution matrix a;

a sound velocity distribution matrix processing unit for aligning the sound velocity distribution matrix c under polar coordinates by using the following formula_sCarrying out normalization and gray level processing;

e (i, j) is for c_s(i, j) normalized and grayed values;

c_s(i, j) is the sound velocity distribution matrix c_sThe value of the element in the ith row and the jth column;

min(c_s) Is the sound velocity distribution matrix c_sMinimum value of (1); max (c)_s) Is the sound velocity distribution matrix c_sMaximum value of (1);

a light absorption energy distribution map gray value conversion unit, which is used for converting the B (i, j) into B '(x, y) in a plane rectangular coordinate according to a formula B' (x, y) ═ B (jcosi, jsini);

a sound velocity profile gray value conversion unit for converting the E (i, j) into E '(x, y) in a plane rectangular coordinate according to the formula E' (x, y) ═ E (jcosi, jsini);

wherein i belongs to [0,360], j belongs to [0, d ];

d is the maximum value of the polar diameter in the theta-l polar coordinate system;

b' (x, y) is a gray value of a point (x, y) in a rectangular coordinate system in the light absorption energy distribution map; e' (x, y) is the grayscale value of the point (x, y) in the rectangular coordinate system in the sound velocity map.

According to the specific embodiment provided by the invention, the invention discloses the following technical effects: the invention obtains the light absorption at the position r of the biological cavityInitial value A of energy receiving amount₀(r) and initial value of sound velocity c_s,0(r), and then calculates the light absorption energy A of the k-th time at the position r, respectively_k(r) and (k + 1) th light absorption energy A_k+1(r), speed of sound c of kth_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r); by judging the absolute difference epsilon of light absorption energy_A,kAnd absolute difference of sound velocityAnd whether the convergence condition is met or not is determined to continue iteration or directly construct a light absorption energy distribution image and a sound velocity distribution image. The invention considers the complexity of biological tissues, comprehensively considers the difference of the speeds, and reconstructs photoacoustic images by combining light absorption energy and different sound velocities, thereby reducing errors caused by uneven sound velocity distribution in tissues to be detected and obtaining high-quality light absorption energy distribution images of the cross section of the cavity.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without creative efforts.

FIG. 1 is a schematic view of a photoacoustic endoscopic imaging method for a biological cavity according to the present invention;

FIG. 2 is a block diagram of a photoacoustic endoscopic imaging system for a biological cavity according to the present invention;

FIG. 3 is a schematic cross-sectional view of a biological cavity in a rectangular XOY coordinate system;

FIG. 4 is a schematic diagram of the multi-layer cavity wall structure in a polar θ -l coordinate system.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

The invention aims to provide a photoacoustic endoscopic imaging method and a photoacoustic endoscopic imaging system for a biological cavity, which are more accurate in constructing sub-images by considering the complexity of biological tissues and combining the conditions of different sound velocities in the process of constructing a light absorption energy distribution image.

In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.

A method of photoacoustic endoscopic imaging of a biological cavity, the method comprising:

obtaining the initial value A of the light absorption energy at the position r of the biological cavity₀(r) and initial value of sound velocity c_s,0(r)。

In this embodiment, the initial value of the light absorption energy is A₀(r) is 0, and the initial value of sound velocity is c_s,0(r)＝1600m/s。

Respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r); k is iteration times, and k is more than or equal to 1;

specifically, this embodiment is based on formula A_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r)) calculating the light absorption energy A of the k-th time_k(r) and the light absorption energy A of the k +1 th order_k+1(r)；

Wherein the content of the first and second substances,

argmin[·]is that g (A)_k(r))+λΦ_TV(A_k(r)) A at the minimum_k(r)；

W(A_k-1(r)) is an approximate blackplug matrix of the first photoacoustic relationship function;

W^-1(A_k-1(r)) is W (A)_k-1(r)) an inverse matrix;

a is a light absorption energy distribution matrix; a is more than or equal to 0, which means the light absorption energy which is more than or equal to 0 in the light absorption energy distribution matrix;

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²is a second photoacoustic relationship function;

p_m(r) is a measurement of the photoacoustic signal at location r;

H(c_s,k-1(r)) is a first operator related to the speed of sound;

p (r) is the theoretical value of the photoacoustic signal at location r;

| | · | is a 2-norm;

Φ_TV(A_k(r)) is a TV regularization term,

λ is the TV regularization parameter;

L_kis about A_k(r) a sparse matrix of variance characteristics;

η > 0 is a constant;

G₁'(r,A_k-1(r),c_s,k-1(r)) is the gradient of the photoacoustic relationship function.

In the specific treatment process, the light energy deposition in the cavity wall tissue under the short pulse laser irradiation is simulated by adopting a Monte Carlo simulation method. Then, solving a discrete photoacoustic wave equation under a polar coordinate system by adopting a finite difference time domain algorithm to obtain a theoretical value of a photoacoustic signal generated by a tissue:

wherein (i, j) is a point on the cross-section of the cavityr coordinates in a theta-l polar coordinate system; Δ θ and Δ l are unit lengths on the θ axis and the l axis, respectively; Δ t is a discrete time interval; n is a discrete time; p is a radical of⁽ⁿ⁺¹⁾(i, j) is a theoretical value of the photoacoustic signal generated by the particle with the position of (i, j) at the time n + 1;andthe vibration speeds of the mass point with the position (i, j) along the theta direction and the l direction at the time n respectively; c. C_s(i, j) is the ultrasonic wave at position (i, j); the propagation velocity of (c); β is the isobaric expansion coefficient of the tissue; c_pIs the specific heat capacity of the tissue; rho₀Is the density of the tissue; i is⁽ⁿ⁾Is the value of the laser pulse function at time n; a (i, j) is the light absorption energy at location (i, j).

Respectively calculating the sound velocity c of the kth time at the position r of the biological cavity_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r)；

In the present embodiment, the first and second electrodes are,

after k iteration is carried out on the third photoacoustic relation function, the sound velocity at the position r is obtained

After the solution is carried out, the solution is obtained,

according to the formulaRespectively calculating the sound velocity c of the k-th time_s,k(r) and the sound speed c of the k +1 st time_s,k+1(r)。

Wherein the content of the first and second substances,

▽f(q_k) Is f (q)_k) A gradient of (a);

f(q_k)＝||p_m(r)-H(q_k)·A_k-1(r)||²；

γ_k-1is the weight parameter, γ, after the k-1 iteration_k-2Is the weight parameter after the k-2 iteration;

weight parameter gamma after the kth iteration_kIs composed of

H(q_k) Is a second operator related to the speed of sound;

d is f (c)_s,k-1(r)) the prestz constant of the derivative;

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²；

c_ssound velocity distribution matrix, c_sAnd > 0 is the sound velocity of 0 or more in the sound velocity distribution matrix.

Calculating the light absorption energy A of the k time_k(r) and the light absorption energy A of the k +1 th order_k+1Absolute difference of light absorption energy ε of (r)_A,kAnd the speed of sound c of the kth time_s,k(r) speed of sound c at the k +1 th time_s,k+1Absolute difference in sound velocity of (r)The convergence condition is epsilon_A,k＜ε_AAnd isWhereinε_AFor the convergence tolerance of the light absorption energy,is the convergence tolerance of the speed of sound.

In this embodiment, the convergence tolerance ε of light absorption energy_ATake 0.01, convergence tolerance of speed of soundTake 0.01.

In the iterative calculation process, if all the iteration values do not meet the convergence condition, the convergence tolerance epsilon of the light absorption energy can be changed_AConvergence tolerance to speed of soundAnd (4) adjusting the convergence condition to complete the iterative process.

Judging the absolute difference epsilon of the light absorption energy_A,kAnd the absolute difference of sound velocityWhether the convergence condition is met or not is judged, and a convergence judgment result is obtained;

if the convergence judgment result is yes, outputting the (k + 1) th light absorption energy A at the position r_k+1(r) and speed of sound c_s,k+1(r) respectively constructing a light absorption energy distribution map and an acoustic velocity distribution map;

if the convergence judgment result is negative, updating the iteration times, and returning to the step of respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r)。

The constructing of the light absorption energy distribution image and the sound velocity distribution image specifically includes:

the light absorption energy distribution matrix A under the polar coordinates is normalized and subjected to gray scale processing by adopting the following formula:

b (i, j) is a value obtained by normalizing and graying a (i, j);

i is the abscissa of the biological cavity position r under the polar coordinate system, and j is the ordinate of the biological cavity position r under the polar coordinate system;

a (i, j) is the element value of the ith row and the jth column in the light absorption energy distribution matrix A;

min (a) is the minimum value of the elements in the light absorption energy distribution matrix a, max (a) is the maximum value of the elements in the light absorption energy distribution matrix a;

adopting the following formula to align the sound velocity distribution matrix c under polar coordinates_sCarrying out normalization and gray level processing;

e (i, j) is for c_s(i, j) normalized and grayed values;

c_s(i, j) is the sound velocity distribution matrix c_sThe value of the element in the ith row and the jth column;

min(c_s) Is the sound velocity distribution matrix c_sMinimum value of (1); max (c)_s) Is the sound velocity distribution matrix c_sMaximum value of (1);

converting the B (i, j) into B '(x, y) in plane rectangular coordinates according to the formula B' (x, y) ═ B (jcosi, jsini);

converting the E (i, j) into E '(x, y) in plane rectangular coordinates according to the formula E' (x, y) ═ E (jcosi, jsini);

wherein i belongs to [0,360], j belongs to [0, d ];

d is the maximum value of the polar diameter in the theta-l polar coordinate system;

b' (x, y) is a gray value of a point (x, y) in a rectangular coordinate system in the light absorption energy distribution map; e' (x, y) is the grayscale value of the point (x, y) in the rectangular coordinate system in the sound velocity map.

The invention also discloses a photoacoustic endoscopic imaging system for the biological cavity, which comprises:

an initial value acquisition module for acquiring an initial value A of light absorption energy at a position r of the biological cavity₀(r) and initial value of sound velocity c_s,0(r)。

In this embodiment, the initial value of the light absorption energy is A₀(r) is 0, and the initial value of sound velocity is c_s,0(r)＝1600m/s。

A light absorption energy calculation module for respectively calculating the light absorption energy A of the kth time at the position r of the biological cavity_k(r) and (k + 1) th light absorption energy A_k+1(r); k is the number of iterations; k is more than or equal to 1.

In this embodiment, the formula A is shown_k(r)＝A_k-1(r)-W^-1(A_k-1(r))G₁'(r,c_s,k-1(r),A_k-1(r)) calculating the light absorption energy A of the k-th time_k(r) and the light absorption energy A of the k +1 th order_k+1(r)。

Wherein the content of the first and second substances,

argmin[·]is that g (A)_k(r))+λΦ_TV(A_k(r)) A at the minimum_k(r)；

W(A_k-1(r)) is an approximate blackplug matrix of the first photoacoustic relationship function;

W^-1(A_k-1(r)) is W (A)_k-1(r)) an inverse matrix;

a is a light absorption energy distribution matrix; a is more than or equal to 0, which means the light absorption energy which is more than or equal to 0 in the light absorption energy distribution matrix;

g(A_k(r))＝||p_m(r)-H(c_s,k-1(r))·A_k(r)||²is a second photoacoustic relationship function;

p_m(r) is a measurement of the photoacoustic signal at location r;

H(c_s,k-1(r)) is a first operator related to the speed of sound;

| | · | is a 2-norm;

Φ_TV(A_k(r)) is a TV regularization term,

λ is the TV regularization parameter;

L_kis about A_k(r) a sparse matrix of variance characteristics;

η > 0 is a constant;

G₁'(r,A_k-1(r),c_s,k-1(r)) is the gradient of the photoacoustic relationship function.

In the specific treatment process, the light energy deposition in the cavity wall tissue under the short pulse laser irradiation is simulated by adopting a Monte Carlo simulation method. Then, solving a discrete photoacoustic wave equation under a polar coordinate system by adopting a finite difference time domain algorithm to obtain a theoretical value of a photoacoustic signal generated by a tissue:

namely, it is

Wherein, (i, j) is the coordinate of a point r on the cross section of the cavity in a theta-l polar coordinate system; Δ θ and Δ l are unit lengths on the θ axis and the l axis, respectively; Δ t is a discrete time interval; n is a discrete time; p is a radical of⁽ⁿ⁺¹⁾(i, j) is a theoretical value of the photoacoustic signal generated by the particle with the position of (i, j) at the time n + 1;andthe vibration speeds of the mass point with the position (i, j) along the theta direction and the l direction at the time n respectively; c. C_s(i, j) is the ultrasonic wave at position (i, j); the propagation velocity of (c); β is the isobaric expansion coefficient of the tissue; c_pIs the specific heat capacity of the tissue; rho₀Is the density of the tissue; i is⁽ⁿ⁾Is the value of the laser pulse function at time n; a (i, j) is the light absorption energy at location (i, j).

A sound velocity calculation module for respectively calculating the positions r of the biological cavitiesThe k-th sound velocity c of_s,k(r) and sound velocity c of (k + 1) th order_s,k+1(r)。

In the present embodiment, the first and second electrodes are,

after k iteration is carried out on the third photoacoustic relation function, the sound velocity at the position r is obtained

After the solution is carried out, the solution is obtained,

according to the formulaRespectively calculating the sound velocity c of the k-th time_s,k(r) and the sound speed c of the k +1 st time_s,k+1(r)。

Wherein the content of the first and second substances,

▽f(q_k) Is f (q)_k) A gradient of (a);

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²；

γ_k-1is the weight parameter, γ, after the k-1 iteration_k-2Is the weight parameter after the k-2 iteration;

weight parameter gamma after the kth iteration_kIs composed of

H(q_k) Is a second operator related to the speed of sound;

d is f (c)_s,k-1(r)) the prestz constant of the derivative;

f(c_s,k-1(r))＝||p_m(r)-H(c_s,k-1(r))·A_k-1(r)||²。

c_ssound velocity distribution matrix, c_sAnd > 0 is the sound velocity of 0 or more in the sound velocity distribution matrix.

An absolute difference calculation module for calculating the light absorption energy A of the kth time_k(r) and the light absorption energy A of the k +1 th order_k+1Absolute difference of light absorption energy ε of (r)_A,kAnd the speed of sound c of the kth time_s,k(r) speed of sound c at the k +1 th time_s,k+1Absolute difference in sound velocity of (r)The convergence condition is epsilon_A,k＜ε_AAnd is

Wherein epsilon_AFor the convergence tolerance of the light absorption energy,is the convergence tolerance of the speed of sound.

In particular, the convergence tolerance ε of the light absorption energy_ATake 0.01, convergence tolerance of speed of soundTake 0.01.

In the iterative calculation process, if all the iteration values do not meet the convergence condition, the convergence tolerance epsilon of the light absorption energy can be changed_AConvergence tolerance to speed of soundAdjusting the convergence condition toThe iterative process is completed.

A convergence judgment result module for judging the absolute difference epsilon of the light absorption energy_A,kAnd the absolute difference of sound velocityWhether the convergence condition is met or not is judged, and a convergence judgment result is obtained;

an image construction module connected with the convergence judgment result module and used for absorbing energy A according to the (k + 1) th light at the position r when the convergence judgment result module is yes_k+1(r) and speed of sound c_s,k+1(r) respectively constructing a light absorption energy distribution map and an acoustic velocity distribution map;

and the iteration updating module is respectively connected with the convergence judgment result module and the light absorption energy calculating module, and is used for updating the iteration times and returning to the light absorption energy calculating module when the convergence judgment result module is negative.

The image construction module specifically comprises:

the light absorption energy distribution matrix processing unit is used for carrying out normalization and gray scale processing on the light absorption energy distribution matrix A under the polar coordinates by adopting the following formula:

b (i, j) is a value obtained by normalizing and graying a (i, j);

i is the abscissa of the biological cavity position r under the polar coordinate system, and j is the ordinate of the biological cavity position r under the polar coordinate system;

a (i, j) is the element value of the ith row and the jth column in the light absorption energy distribution matrix A;

min (a) is the minimum value of the elements in the light absorption energy distribution matrix a, max (a) is the maximum value of the elements in the light absorption energy distribution matrix a;

a sound velocity distribution matrix processing unit for aligning the sound velocity distribution matrix c under polar coordinates by using the following formula_sPerforming normalization and gray scale processingC, processing;

e (i, j) is for c_s(i, j) normalized and grayed values;

c_s(i, j) is the sound velocity distribution matrix c_sThe value of the element in the ith row and the jth column;

min(c_s) Is the sound velocity distribution matrix c_sMinimum value of (1); max (c)_s) Is the sound velocity distribution matrix c_sMaximum value of (1);

a light absorption energy distribution map gray value conversion unit, which is used for converting the B (i, j) into B '(x, y) in a plane rectangular coordinate according to a formula B' (x, y) ═ B (jcosi, jsini);

a sound velocity distribution image gradation value conversion unit for converting the E (i, j) into E '(x, y) in plane rectangular coordinates according to a formula E' (x, y) ═ E (jcosi, jsini);

wherein i belongs to [0,360], j belongs to [0, d ];

d is the maximum value of the polar diameter in the theta-l polar coordinate system;

b' (x, y) is a gray value of a point (x, y) in a rectangular coordinate system in the light absorption energy distribution map; e' (x, y) is the grayscale value of the point (x, y) in the rectangular coordinate system in the sound velocity map.

For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.

The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.