JP6799321B2 - Optical ultrasonic imaging device and method, control program of optical ultrasonic imaging device, and recording medium - Google Patents

Description

The present invention relates to an optical ultrasonic imaging apparatus and method for imaging an object to be inspected by using an optical ultrasonic imaging method, a control program of the optical ultrasonic imaging apparatus, and a recording medium.

Photoacoustic imaging (hereinafter also referred to as PA imaging or PA imaging) is a promising new technique for early clinical diagnosis of cancer, tumor vascular neoplasms, and many other diseases (eg, non-PA imaging). See Patent Document 1). PA imaging is useful in many medical applications because it can non-invasively visualize blood vessels in a living body with high spatial resolution without a contrast medium. Here, PA imaging utilizes the photoacoustic effect. That is, the imaged tissue (ie, blood vessels) absorbs the laser energy of short pulse near-infrared irradiation and converts it into heat, resulting in thermoelastic expansion and emission of ultrasonic waves. The three-dimensional structure of an object can be reconstructed by detecting the emitted photoacoustic wave (see, for example, Non-Patent Document 2).

The reconstruction problem has the following two major difficulties. As a first difficulty, the data acquisition speed of PA imaging is limited by the laser repetition rate and the number of ultrasonic sensors. The number of samplings is limited for real-time imaging, reducing system costs. This inadequate sampling usually causes streaking artifacts. As a second difficulty, assuming that the speed of sound is correctly given in most PA reconstruction methods, the accuracy of reconstruction largely depends on the accuracy of the given velocity. The speed of sound often changes with changes in temperature and tissue, and it is essential to adjust accordingly.

For example, in Patent Document 1, a subject for acquiring highly accurate subject information by using a probe configured by arranging a large number of oscillators along a substantially spherical cap shape or a substantially spherical band shape. The information acquisition device is disclosed. This subject information acquisition device includes a plurality of oscillators that receive acoustic waves from a measurement target and generate an electric signal, and a support that supports so that at least a part of the directional axes of the plurality of oscillators are gathered. From a point sound source, a means for controlling the relative position between the measurement target and the support, a processing means for obtaining the distance from each of a plurality of oscillators to the attention position of the measurement target and generating characteristic information at the attention position, and a point sound source It is characterized by using a subject information acquisition device having a correction means for obtaining correction data of position information of each of the plurality of oscillators by using the reception time of an acoustic wave to the plurality of oscillators and the sound velocity.

Further, Patent Document 2 discloses that an ultrasonic imaging apparatus capable of acquiring an ultrasonic image having a higher resolution than the conventional one is provided by optimizing the sound velocity used for ultrasonic imaging. .. In the ultrasonic imaging apparatus, the resolution optimization unit determines the optimum sound velocity according to the tissue component for each position in the scanning cross section, and uses this optimum sound velocity for each received beam from each position in the scanning cross section. The reception delay time of is calculated. The control processor uses the reception delay time calculated using the optimum speed of sound to perform delay addition processing in scanning to acquire an ultrasonic image actually used for diagnosis. As a result, it is possible to correct the deviation between the set sound velocity used for calculating the reception delay time and the actual in-vivo sound velocity, and acquire an ultrasonic image with optimized resolution.

Japanese Unexamined Patent Publication No. 2016-055162 Japanese Unexamined Patent Publication No. 2008-264531

Kitai, T. et al., "Photoacoustic mammography: initial clinical results," Breast Cancer, 21, pp.146-153 (2014). Li C. et al. "Photoacoustic tomography and sensing in biomedicine," Physics in Medicine and Biology, 54 (19), pp.59-97 (2009). Candes, E.J et al., "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information," IEEE Transactions on Information Theory, pp. 489-509, (2006). Zhang, Z. et al., "Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction," Ultrasonics, 52, pp.1046-1055, (2012). Xu M et al., "Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries, IEEE Trans on Biomedical Engineering, 50 (9), pp.1086-1099, (2003). Xu Y et al., "Exact frequency-domain reconstruction for thermoacoustic tomography I: Planar geometry," IEEE Medical Imaging, 21 (7), pp.823-828, (2002). Xu, M. et al., "Universal back-projection algorithm for photoacoustic computed tomography," Physical Review, 71, 016706, (2005). Guo, Z. et al., "Compressed sensing in photoacoustic tomography in vivo," Journal of Biomedical Optics, 15 (2), 021311, (2010). Treeby, BE. Et al., "Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach," Journal of Biomedical Optics, 16 (9), 090501, (2011). Li, C. et al., "An Efficient Augmented Lagrangeian Method with Applications to Total Variation Minimization," Computational Optimization and Applications, 56 (3), pp.507-530, (2013). Treeby, BE. Et al., "K-wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields," Journal of Biomedical Optics, 15 (2), 021314, (2010). Kruger, RA. Et al., "Dedicated 3D photoacoustic breast imaging," Medical Physics, 40 (11), 113301, pp.1-8 (2013). Mast, TD. Et al., "Simulation of ultrasonic pulse propagation through the abdominal wall," Journal of Acoustic Society of America, 102 (2), pp.1177-1190, (1997). Kazufumi Mimura, "Compressed Sensing-Reconstruction of Sparse Information and Its Algorithms", Research Institute for Mathematical Sciences, Vol. 1803, 2012 26-56, (2012). R. Baraniuk et al., “Random projections of smooth manifolds,” 11 Foundations of Computational Mathematics, Vol. 9, No. 1, pp. 51-77, 2009. Y. C. Eldar et al., "Robust recovery of signals from a union of subspaces," IEEE Transactions on Information Theory, Vol. 55, No. 11, pp. 5302-5316, 2009. R. Baraniuk, V. et al., "Model-based compressive sensing," Preprint, 2008 A. W. Fitzgibbon, "Robust Registration of 2D and 3D Point Sets," Proceedings of the British Machine Vision Conference, 2001. Rueckert, D. et al., "Nonrigid registration using free-form deformations: application to breast MR images," IEEE Transactions on Medical Imaging, 8, pp.712-721 (1999). R. Bise et al., "Vascular Registration in Photoacoustic Imaging by Low-Rank Alignment via Foreground, Background and Complement Decomposition," Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention, pp.326-334, (2016) ..

However, when the device of Patent Document 1 or 2 is used to image an inspection object, the speed of sound needs to be known, but the speed of sound fluctuates depending on the structure or temperature of the inspection object, and is assumed. Unlike the speed of sound, there was a problem that the image deteriorated and the accuracy of the image was still low.

An object of the present invention is to solve the above problems and to provide an optical ultrasonic imaging device and method capable of imaging an object with higher accuracy than the prior art, a control program of the optical ultrasonic imaging device, and a recording medium. To provide.

The optical ultrasonic imaging apparatus according to the first invention irradiates an object with light by an optical ultrasonic imaging method, and uses a plurality of sensors to generate ultrasonic waves from the object as a plurality of detection signals. In an optical ultrasonic imaging apparatus provided with a control means for performing an imaging process for imaging the object based on sound pressure data including a plurality of detection signals indicating the sound pressure spatial distribution after detection.
The control means is an objective function relating to the initial sound pressure vector when the multiplication result obtained by multiplying the projection matrix of the sound pressure data by the initial sound pressure vector is used as the sound pressure data based on the predetermined initial sound pressure data. The value of the objective function representing the sparseness of the image data of the object or the sparseness of the base data of the conversion when a predetermined conversion is performed on the image data is a predetermined threshold value. The imaged image data based on the sound velocity vector is calculated by simultaneously estimating the sound pressure vector and the sound velocity vector by changing the sound velocity vector including the sound velocity in a plurality of regions having different sound velocity as follows. It is characterized by that.

In the optical ultrasonic imaging device, the control means calculates imaged image data based on the sound velocity vector so that the value of the objective function is substantially minimized.

Further, in the optical ultrasonic imaging apparatus, the objective function is the L1 norm, L2 norm, of the image data of the object or the data of the basis of the conversion when a predetermined conversion is performed on the image data. Alternatively, it is characterized by total variation or a combination thereof.

Further, in the optical ultrasonic imaging device, the control means outputs the image data of the object in the process of simultaneously estimating the sound pressure vector and the sound velocity vector, or the image data after the imaging to the display unit. It is characterized by being visualized by.

Furthermore, in the optical ultrasonic imaging device, the control means changes at least one of the spatial resolution and the temporal resolution of the optical ultrasonic imaging method in ascending order, and uses the optical ultrasonic imaging method. Each of the imaging processes is executed, and the sound velocity vector obtained by each imaging process is used as the next initial sound velocity vector to perform the imaging process.

Furthermore, in the optical ultrasonic imaging device,
The object includes first and second regions having different sound velocities from each other.
The control means
(1) Using a predetermined initial sound velocity vector, the object is imaged using a predetermined imaging method, and image data is calculated.
(2) The second region is estimated by performing predetermined image processing on the calculated image data.
(3) The initial sound velocity map is calculated by aligning the estimated second region with the predetermined sound velocity initial value atlas showing the initial sound velocity value of each region of the object by the predetermined alignment method. And
(4) The imaging process according to any one of claims 1 to 5 is performed using the aligned initial sound velocity map.

The optical ultrasonic imaging method according to the second invention is an optical ultrasonic imaging method in which an object is irradiated with light and ultrasonic waves generated from the object are used as a plurality of detection signals by a plurality of sensors. Optical ultrasound for an optical ultrasound imaging device equipped with a control means that performs imaging processing to image the object based on sound pressure data consisting of a plurality of detection signals indicating the sound pressure spatial distribution. In the imaging method
When the control means obtains the sound pressure data as the result of multiplying the projection matrix of the sound velocity data by the initial sound pressure vector based on the predetermined initial sound pressure data, the objective function relating to the initial sound pressure vector. The value of the objective function representing the sparseness of the image data of the object or the sparseness of the base data of the conversion when a predetermined conversion is performed on the image data is a predetermined threshold value. Imaged image data based on the sound velocity vector is calculated by simultaneously estimating the sound pressure vector and the sound velocity vector by changing the sound velocity vector including the sound velocity in a plurality of regions having different sound velocity as follows. It is characterized by including steps.

In the optical ultrasound imaging method, the control means includes a step of calculating imaged image data based on the sound velocity vector so that the value of the objective function is substantially minimized. To do.

Further, in the optical ultrasonic imaging method, the objective function is the L1 norm, L2 norm, of the image data of the object or the data of the basis of the conversion when a predetermined conversion is performed on the image data. Alternatively, it is characterized by total variation or a combination thereof.

Further, in the optical ultrasonic imaging method, the control means outputs the image data of the object in the process of simultaneously estimating the sound pressure vector and the sound velocity vector, or the image data after the imaging to the display unit. It is characterized by including a step of visualizing by.

Furthermore, in the optical ultrasonic imaging method, the control means changes at least one of the spatial resolution and the temporal resolution of the optical ultrasonic imaging method in ascending order, and the optical ultrasonic imaging method is performed. Each of the imaging processes is executed, and the sound velocity vector obtained by each imaging process is used as the next initial sound velocity vector to perform the imaging process.

Furthermore, in the optical ultrasonic imaging method,
The object includes first and second regions having different sound velocities from each other.
The control means
(1) Using a predetermined initial sound velocity vector, the object is imaged using a predetermined imaging method, and image data is calculated.
(2) The second region is estimated by performing predetermined image processing on the calculated image data.
(3) The initial sound velocity map is calculated by aligning the estimated second region with the predetermined sound velocity initial value atlas showing the initial sound velocity value of each region of the object by the predetermined alignment method. And
(4) It is characterized in that the step of performing the imaging process according to any one of claims 1 to 5 is included using the aligned initial sound velocity map.

A computer-executable control program for an optical ultrasonic imaging device according to a third aspect of the invention is characterized by including each step of the optical ultrasonic imaging method.

A computer-readable recording medium according to a fourth aspect of the invention is characterized by storing a control program of the optical ultrasonic imaging apparatus.

Therefore, according to the optical ultrasonic imaging apparatus or the like according to the present invention, it is possible to image an inspection object with higher accuracy than in the prior art.

It is a schematic cross-sectional view and block diagram which shows the structural example of the optical ultrasonic imaging system which concerns on one Embodiment of this invention. The result image of the optical ultrasonic imaging process, (a) is a true value image, (b) is an image by the back projection method, and (c) is a compressed sensing method using an inaccurate sound velocity (hereinafter, CS method). The image according to () and (d) are the images by the CS method using an accurate sound velocity. It is a top view which shows the ultrasonic wave propagation at time t1 in the case of a constant sound velocity c. It is a flowchart which shows the optical ultrasonic imaging process which concerns on Embodiment 1. It is a front view which shows the relationship between the true value image and the position of a plurality of sensors 22. It is a graph which shows the PSNR (Peak Signal to Noise Ratio) characteristic with respect to the sound velocity by the back projection method and CS method which are conventional examples, and the imaging method which concerns on Embodiment 1. It is a graph which shows the estimated sound velocity with respect to the sound velocity of the data set in the true value image and the image by the imaging method which concerns on Embodiment 1. FIG. Images at various sound velocities, (a) is an image by the back projection method, (b) is an image by the CS method, and (c) is an image by the imaging method according to the first embodiment. .. (A) is an imaging result image using the combination of sound velocities, (b) is an image by the back projection method, (c) is an image by the CS method, and (d) is an image by imaging according to the first embodiment. , (E) is an enlarged image of the image of (c), and (f) is a diagram showing an enlarged image of the image of (d). It is a flowchart which shows an example of the high precision optical ultrasonic imaging process which concerns on Embodiment 2 of this invention. It is a flowchart which shows an example of the high precision optical ultrasonic imaging process which concerns on Embodiment 3 of this invention. It is a flowchart which shows the imaging process using the sound velocity initial value atlas which concerns on Embodiment 4 of this invention. It is a schematic diagram which shows the sound velocity initial value atlas used in the imaging process of FIG. It is a schematic diagram which shows the process in step S31 of FIG. It is a schematic diagram which shows the process in step S32 of FIG. It is a schematic diagram which shows the process before the alignment of step S33 of FIG. It is a schematic diagram which shows the process after the alignment of step S33 of FIG. It is a schematic diagram which shows the process in step S34 of FIG. It is a schematic diagram which shows the process in step S35 of FIG.

Hereinafter, comparative examples and embodiments according to the present invention will be described with reference to the drawings. In each of the following embodiments, the same reference numerals are given to the same components.

An object of the embodiment of the present invention is a new compression for imaging while estimating the speed of sound at the same time when recovering the three-dimensional structure of an imaged inspection object (object) from a limited number of sampled data. The purpose is to propose an imaging method by a sensing method (hereinafter referred to as a CS method). The CS method is based on the basic fact that many signals are represented by a suitable reference or dictionary with only a few non-zero coefficients, to improve reconstruction accuracy even from a limited number of observations. It is valid. Here, nonlinear optimization makes it possible to recover such a signal from a small number of measurements (see, eg, Non-Patent Document 3). Imaging by the CS method has been applied to various medical imaging methods, but there are various methods such as MRI, CT, PA and the like (see, for example, Non-Patent Document 4). The current CS method of PA imaging has a problem that it does not function effectively in the case of an incorrect sound velocity.

FIG. 2 is an image obtained as a result of optical ultrasonic imaging processing, FIG. 2 (a) is a true value image, FIG. 2 (b) is an image obtained by the back projection method, and FIG. 2 (c) uses an inaccurate sound velocity. The image obtained by the CS method and FIG. 2 (d) are diagrams showing an image obtained by the CS method using an accurate sound velocity.

As is clear from FIG. 2D, when the correct speed of sound is given, the CS method can obtain the same good results as the true value image (FIG. 2A). In contrast, if the speed of sound given is incorrect, both the back projection method and the CS method will not give very good images. (See FIGS. 2 (b) and 2 (c)).

In order to solve the above problems, in the present embodiment, the reconstruction problem can be simultaneously solved from the ultrasonic signal having sparsity, and the image processing by the new CS method for PA imaging capable of estimating the sound velocity. To propose. Various CS method optimization techniques can be used within the embodiment. In particular, we used a method of optimizing the CS method that minimizes total variation using the extended Lagrange multiplier method. The results of the experiment by the inventor were able to demonstrate the effectiveness of this proposed method, evaluate the speed of sound, and significantly improve the image quality, and the present embodiment is based on the speed of sound. This is the first attempt to tackle the problem of image reconstruction while simultaneously estimating. In the conventional example, this problem has not been solved even assuming a uniform sound velocity.

FIG. 1 is a schematic cross-sectional view and a block diagram showing a configuration example of an optical ultrasonic imaging system according to an embodiment of the present invention.

In FIG. 1, the optical ultrasonic imaging system according to the present embodiment includes an optical ultrasonic light source, a detection device 20, and an optical ultrasonic imaging device 10.

In the optical ultrasonic light source and the detection device 20, for example, a plurality of sensors 22 for detecting ultrasonic waves are provided on a hemispherical sensor support portion 21 at a predetermined distance from each other. Water 24 is immersed in the hemispherical surface to support the mesh plate 23 on which the inspection object 1 is placed. It is assumed that the inspection object 1 is submerged in water 24. A laser light source 25 for radiating a laser beam having a predetermined wavelength for PA imaging to the inspection object 1 is provided at the bottom of the sensor support portion 21. Here, the laser light source 25 is controlled by a controller 11 described later, and the ultrasonic waves radiated from the inspection object 1 in response to the laser light are detected by the plurality of sensors 22 by PA imaging. The ultrasonic waves detected by the plurality of sensors 22 are converted into electrical signals and then sent to a signal receiving circuit 12 described later.

The inspection object 1 is an imaging object. Specific examples include a living body such as a breast, and a phantom that simulates the acoustic and optical characteristics of a living body in adjusting a device. The acoustic characteristics are specifically the propagation velocity and the attenuation rate of the acoustic waves, and the optical characteristics are specifically the light absorption coefficient and the scattering coefficient. It is preferable that the mesh plate 23 supports the inspection object 1 and has transmission characteristics for light and ultrasonic waves.

The optical ultrasonic imaging device 10 includes a controller 11 that controls the operation of the device, a signal receiving circuit 12, a data memory 13, an operation unit 14, a display unit 15, and an image processing unit 16. To. The signal receiving circuit 12 receives a plurality of electric signals converted by the plurality of sensors 22, A / D-converts them at a predetermined sampling frequency, and stores them in the data memory 13. The operation unit 14 includes a mouse and a keyboard for the operator of the optical ultrasonic imaging device 10 to input data or instructions. The display unit 15 displays the input data or instructions, and visualizes and displays the result image during or after the optical ultrasonic imaging process. The image processing unit 16 is controlled by the controller 11 and executes optical ultrasonic imaging processing described in detail later.

Although the hemispherical sensor support portion 21 is used in the embodiment of FIG. 1, the present invention is not limited to this, and the directional axes of the plurality of sensors 22 are gathered in a predetermined region to form a predetermined detection region. Any arrangement can be used.

In the photoultrasound imaging system of FIG. 1, photoacoustic effects can be used to non-invasively visualize, for example, living blood vessels in cancer and many others. It is useful for early clinical diagnosis of the disease. The tissue to be imaged, which is the object to be inspected 1, absorbs laser energy, converts the energy into thermoelastic expansion, and emits ultrasonic waves. By detecting optical ultrasonic waves, the three-dimensional structure of the object of the inspection object 1 can be reconstructed. In most PA imaging, the speed of sound is assumed to be given correctly, but it is known that the speed of sound varies depending on the temperature of the propagation medium and tissue. Therefore, it is difficult to know this value accurately in a living body, but in the present embodiment, a new compressed sensing method for PA that simultaneously solves the reconstruction problem and estimates the speed of sound is proposed.

For example, in the literature related to this embodiment, the PA reconstruction method based on the inverted spherical Radon transform is proposed in both the time domain (see, eg, Non-Patent Document 5) and the frequency domain (see, eg, Non-Patent Document 6). Has been done. The filtered back projection method proposed in Non-Patent Document 7 is the most popular because of its convenience. One problem with these methods is the serious artifacts that are generated when the number of samples is inadequate. In order to improve the image quality using sampling data having sparsity, the CS theory was applied to the reconstruction by the CS method. Non-Patent Document 8 demonstrated the possibility of using CS theory in PA reconstruction. In addition, Non-Patent Document 4 proposes a method of iteratively minimizing total variation (TV) as an objective function due to gradient descent. These methods assume that the speed of sound is given correctly. To select the speed of sound, Non-Patent Document 9 proposes an autofocus approach that repeatedly reconstructs an image using various candidate sound velocities, but if the speed of sound is assumed to be spatially uniform,
Then select the best speed of sound that maximizes the sharpness metric. However, when there are a plurality of regions having different sound velocities, there is a problem that the number of combinations increases exponentially. As far as the inventors know, this challenging problem is addressed first, and the method according to this embodiment is characterized by iteratively updating the set of sound velocities in optimization.

First, compressed sensing for PA imaging by sound velocity estimation will be described below. Here, we show how the initial sound pressure, which briefly describes the discretized formulation of PA propagation, propagates to the detector as ultrasonic waves. In PA imaging, the reconstructed target image map can be thought of as the spatial distribution of acoustic initial sound pressure.

FIG. 3 is a plan view showing ultrasonic wave propagation at time t1 in the case of a constant sound velocity c. In FIG. 3, the ultrasonic signals from the light absorbing objects A and C have just arrived at another sensor 22 at time t1, and the distance R can be calculated from the product of the speed of sound c and the time t1. The observation signal of the sensor 22 at time t1 can be considered to be the total of the ultrasonic signals that have just arrived. According to Non-Patent Document 4, this forward propagation problem is expressed by the following equation.

b = A (c) u (1)

Here, TK × 1 vector b represents an electric signal measured by a sensor array including a plurality of sensors 22, K is the number of sensors 22, and T is the number of time slots. The vector u represents the spatial distribution of the initial sound pressure, where the magnitude of u (N = N _x × N _y × N _z ) depends on the visual field. N _x is the number of sensors 22 in the x direction, N _y is the number of sensors 22 in the y direction, and N _z is the number of sensors 22 in the z direction.

Here, the TK × N matrix A (c) is defined as a projection matrix, and how the sound pressure propagates to each position and time is described. Examples of A (c) are described in Supplementary Explanation A2.

The reconstruction problem includes solving the distribution of the initial sound pressure u satisfying b = A (c) u. If the sampled data is inadequate, the opposite problem is a bad idea. This is because it results in striped artifacts in the reconstructed image. To address this issue, the CS method takes into account the characteristics of the target image in optimizing the reconstruction problem. Here, re-in the CS method, assuming sparseness in a more basic domain, such as a Fourier domain (see, eg, Non-Patent Document 3), small total variation (see, eg, Non-Patent Document 4). The configuration problem can be formulated as follows.

(2)

Here, "s.t." is an abbreviation for "so that" and indicates a constraint condition. f is an objective function that describes the characteristics (features such as sparsity) of the target image. Typically, the L1 norm | Ψu | ₁ or total variation TV (u) is used for the objective function f. Where Ψ is a sparse transformation matrix. If the sound velocity map c is given correctly, the projection matrix A can be calculated and fixed. In the above equation (2), a sound such that the value of the objective function f (u) becomes smaller by changing the spatial distribution u of the sound pressure and the speed of sound c under the condition of b = A (c) u. The spatial distribution u of pressure and the speed of sound c are obtained. Therefore, the image can be reconstructed by using the CS method (see, for example, Non-Patent Document 8). However, solving this problem is difficult to estimate u and c at the same time because there are many unknowns and the projection matrix A (c) is non-linear at the normal speed of sound c.

Next, iterative linearization for estimating the speed of sound c will be described below.

As mentioned above, the optimization problem of equation (2) is difficult to handle directly. Here, one of the major difficulties is the large number of 2N unknown parameters. To address this problem, the inventors can mitigate it by introducing additional but realistic assumptions. In the imaging target, the present inventors have only a few regions having different sound velocities, and the positions of these regions are generally known. That is, for example, since the position of the tissue of the test object 1 does not differ greatly from patient to patient, a small difference in position does not have a serious effect on the reconstruction. Here, the sound velocity vector c is expressed by the following equation.

c = [c ₁ , c ₂ , ..., c _L ] (3)

Here, cl is the speed of sound in the l-th region, and L is the number of regions having different velocities. Hereinafter, the present inventors estimate L velocities (sound velocities) in the reconstruction problem. Another difficulty is non-linearity, where the acoustic propagation function A (c) is not linear with respect to the speed of sound c. When the change in sound velocity c is small, the projection matrix A (c) can be approximated to Δc = [Δc ₁ | ... | c _L ] by linearizing the current estimated value of sound velocity c. it can. Here, the elements of A _{i, j, k} (c + Δc) can be approximated as in the following equation.

(4)

This leads to the optimization problem in the unknowns u and Δc, and the objective function can be described as follows.

(5)

Since linearization is only retained locally, one should not expect the solution (c + Δc) from equation (5) to accurately solve equation (1). In order to find the local minimum value of equation (1), we can iteratively linearize the current estimate of the sound velocity vector c and solve a series of optimizations of equation (5). As can be seen from the experimental results described later, this repetition can estimate the initial sound pressure u and the speed of sound c as long as the speed difference of the initial sound velocity is not too large. The complete optimization procedure is summarized as the optical ultrasound imaging process of FIG. In the equation (5), the result data (u, Δc) when the value of the objective function is the minimum value is obtained, but the present invention is not limited to this, and is substantially the minimum value or the minimum value. The above-mentioned result data which is close to and is a value equal to or less than a predetermined threshold value may be obtained.

FIG. 4 is a flowchart showing the optical ultrasonic imaging process according to the first embodiment. The optical ultrasound imaging process represents a linearized sequence of optimizations.

In step S1 of FIG. 4, first, the initial sound velocity vector is c. Next, in step S2, let b be the vectorized detection signal data.

Here, b is vectorized detection signal data in a space formed by T time slots and K sensors.

In step S3, the projection matrix A (c) and the Jacobian matrix are calculated using the following equations based on the current sound velocity vector c.

(6)

Next, in step S4, the result data (u ^* , Δc ^* ) is calculated by performing the calculation of the linear optimization using the following equation. Here, Δc is a sound velocity change amount vector, and Δc ^* is a sound velocity change amount vector after linear optimization.

(7)

In the formula (7), specifically, the formula (10) in the supplementary explanation A1 is used.

In step S5, c + Δc ^* is calculated and set as the next sound velocity vector c. In step S6, the signal data b = A (c) u ^* is calculated using the current sound pressure data u ^* , the signal data b is calculated and output to the display unit 15, and the intermediate imaged image is displayed on the display unit. Display on 15. Then, in step S7, it is determined whether convergence conditions (Δc ^* <Δc _th) satisfies, while proceeds step S8 when YES, the if NO the process returns to the step S3. Here, .DELTA.c _th is the threshold of sound change amount vector. Further, in step S8, the result data (u ^* , Δc ^* ) and the signal data b are output to the display unit 15, the imaged image is displayed, and the process is completed.

Next, the optimization of each iteration will be described below.

In the proposed optimization method, various CS methods can be used. In this embodiment, the objective function of total variation (TV) is used. The TV minimization method, which minimizes the objective function of total variation, is widely used in modern biomedical image reconstruction and is also applied to PA imaging (see, for example, Non-Patent Document 4). The objective function of total variation is premised on the sparse gradient of the underlying signal (having sparsity). This is to make the reconstructed image sharper by preserving the edges or boundaries more accurately compared to the minimization of l1. In order to solve the minimization problem using the objective function of total variation, the present inventors use a method (extended Lagrange multiplier method) including an extended Lagrange multiplier (ALM) (see, for example, Non-Patent Document 10). Specifically, the discrete TV model for CS reconstruction is expressed by the following equation.

(8)

here
The sum of indicates a period in the objective function of total variation, D i _u is a discrete gradient of the sound pressure space distribution u in the voxel i. The optimization problem of equation (8) can be solved by minimizing the extended Lagrange function of equation (8).

(9)

Here, D is N × N matrix, i-th row is D _i, the Σ and Δ is the Lagrange multiplier matrix for each constraint, the β and μ is a positive scalar. For details on optimization, refer to Supplementary Explanation A1.

Next, the experimental evaluation will be described below.

The present inventors performed a simulation to verify the effectiveness of this method for PA reconstruction. Photoacoustic propagation was simulated using MATLAB's k-wave toolbox (see, eg, Non-Patent Document 11). This toolbox is widely used for photoacoustic simulation. Here, the simulation was performed in two dimensions.

FIG. 5 is a front view showing the relationship between the true value image and the positions of the plurality of sensors 22. In FIG. 5, the gray scale at the right end indicates the strength of the initial acoustic sound pressure. The black dots indicate the positions of the sensor 22. Here, the number of sensors 22 is 70. The size of the reconstructed image was 216 × 216 pixels, and the size of each pixel was 0.0463 [mm].

First, we evaluated the method when the speed of sound is constant. Because test datasets were generated using various sound velocities from 1460 to 1550 [m / s], the speed of sound increment was 10 [m / s], and the speed of sound in living tissues was often within this range. Is. The results are obtained by a back projection method using a time reversal method (see, for example, Non-Patent Document 11) and a CS method (see, for example, Non-Patent Document 10) that minimizes the total variation objective function using ALM. Compared with the results. These methods were applied to each simulated dataset with a speed of sound of 1500 [m / s].
The speed in the reconstruction process was estimated with 1500 [m / s] set as the initial speed of sound.

8 (a) is an image obtained by a back projection method, FIG. 8 (b) is an image obtained by a CS method, and FIG. 8 (c) is an image according to the first embodiment. It is a figure which shows the image by the chemical conversion method. Here, it can be seen that the back projection and CS methods suffer artifacts due to the speed of sound given incorrectly. In contrast, the method according to this embodiment was able to successfully reproduce an image in the sound velocity range of 1480 to 1540 [m / s]. In order to quantitatively evaluate the image quality, the peak signal-to-noise ratio (PSNR) of the reconstructed image was calculated with reference to the original phantom. Here, the greater the value of PSNR, the better the quality of the image.

FIG. 6 is a graph showing PSNR (Peak Signal to Noise Ratio) characteristics with respect to sound velocity by the back projection method and the CS method, which are conventional examples, and the imaging method according to the first embodiment. As is clear from FIG. 6, the image quality of the result according to the present embodiment is the highest among 1480 to 1540 [m / s].

FIG. 7 is a graph showing the estimated sound velocity with respect to the sound velocity of the data set in the true value image and the image obtained by the imaging method according to the first embodiment. That is, FIG. 7 is an evaluation of the sound velocity estimation of the present embodiment. The method according to this embodiment normally estimated an accurate speed of sound from 1480 to 1540 [m / s]. Next, we evaluated the method when there are multiple regions with different velocities.

9 (a) is an imaging result image using the combination of sound velocities, FIG. 9 (b) is an image by the back projection method, FIG. 9 (c) is an image by the CS method, and FIG. 9 (d) is an image. 9 (e) is an enlarged image of the image of FIG. 9 (c), and FIG. 9 (f) is an enlarged image of the image of FIG. 9 (d) according to the first embodiment. Here, a test data set was generated using the combination of sound velocities shown in FIG. 9 (a), and the range of sound velocities was 1470 to 1540 [m / s]. As is clear from FIG. 9, artifacts occur in the results of the back projection method (FIG. 9 (b)) and the CS method (FIG. 9 (c)) when the speed of sound is 1540 to 1510 [m / s]. In contrast, the method according to this embodiment restored images normally and estimated the speed of sound for all combinations. Here, the average PSNR was as follows.
(1) Back projection method: 21.75;
(2) CS method: 23.46;
(3) Method according to the first embodiment: 28.45,
The average sound estimation error was 2.17 [m / s] (1.4%).

As described above, according to the method according to the present embodiment, a new reconstruction method for PA imaging is proposed, the reconstruction problem is solved from the ultrasonic signal having sparsity, and the sound velocity is simultaneously controlled by using the CS method. I was able to solve it by estimating. The results of the experiment demonstrated the effectiveness of the method according to the embodiment and succeeded in significantly improving the speed of sound and image quality.

(Supplemental Information A1) Optimization of ALM.
The optimization problem of equation (9) is iteratively minimized as in equation (10) below.

(10)

Here, "sgn" is a sign function, k is an index of iterations, and is determined using a non-monotone line search such as Amijo.

(Supplementary explanation A2) Projection matrix.
According to Non-Patent Document 13, an attenuated sine wave of the following equation is used for the projection matrix A (c).

(11)

Here, ω ₀ is the frequency of the ultrasonic signal, i is the position index, t _j is the time of the jth time slot, and _dl is the speed of sound in the path between the kth sensor 22 and the sound source. is a distance on the areas with c _l.

The following may be used for the projection matrix A (c).
Average
, The time slot t _j of the Gauss distribution function, etc. with the variance σ
A function with a shape that maximizes when

(Modified Example of Embodiment 1)
As the objective function used in the above first embodiment, the following various objective functions, the objective function, or a combination thereof can be used.
(1) When the sparsity in the base domain (base domain) after the sound pressure spatial distribution u representing the reconstructed image is converted by a predetermined conversion function is used:
L1 norm | Ψu | ₁ is used. Here, for example, the basis of the following transformation function can be used.
(A) Fourier transform: Ψ is the basis of the inverse Fourier transform.
(B) Wavelet transform: Ψ is the basis of the inverse Wavelet transform.
(C) Discrete cosine transform (d) Dictionary of orthogonal bases (see, for example, Non-Patent Document 16).
(2) When using the sparsity of the sound pressure spatial distribution u itself representing the reconstructed image:
The L1 norm | u | ₁ is used. Here, for example, the following function is used.
(A) TV (total variation): total variation;
(B) L1 norm: | TV (u) | ₁ : Local change in brightness is sparse;
(C) L2 norm: | TV (u) | ₂ : The total amount of changes in local brightness is small;
(D) low dimensional smartphone manifolds (see, for example, Non-Patent Document 17);
(E) union of sub-spaces (see, for example, Non-Patent Document 18);
(F) block-sparseness (see, for example, Non-Patent Document 19).
(3) A combination of the above (1) and (2) (for example, a linear combination) may be used.

Embodiment 2.
FIG. 10 is a flowchart showing an example of high-precision optical ultrasonic imaging processing according to the second embodiment of the present invention. In the second embodiment, the spatial resolution is increased in ascending order by changing the number N of the sensors 22 from a small number to a large number, for example, N = 20 × 20 × 20, N = 40 × 40 × 40, N = 100 × 100 × 100. The optical ultrasonic imaging process (S100) according to the first embodiment is executed, and the sound velocity vector c ^* obtained in each process is used as the next initial sound velocity vector c to perform the imaging process. It is characterized by that.

In FIG. 10, in step S11, the number of sensors N = 20 × 20 × 20, and the optical ultrasonic imaging process of step S100 is executed in step S12 to obtain the sound velocity vector c ^* . Next, in step S13, the sound velocity vector c ^* obtained above is used as the initial sound velocity vector for the next imaging process.

Next, in step S14, with the number of sensors N = 40 × 40 × 40, the optical ultrasonic imaging process of step S100 is executed in step S15 to obtain the sound velocity vector c ^* . Next, in step S16, the sound velocity vector c ^* obtained above is used as the initial sound velocity vector for the next imaging process.

Next, in step S15, assuming that the number of sensors N = 100 × 100 × 100, the optical ultrasonic imaging process of step S100 is executed in step S16, and the result data (u, c _* ) including the sound velocity vector c ^* is obtained. Is output to the display unit 15 for visualization and display.

As described above, according to the second embodiment, the spatial resolution is changed in ascending order, the optical ultrasonic imaging process (S100) according to the first embodiment is executed, and the sound velocity vector c obtained in each process is performed. The accuracy of the imaged image can be improved by performing the imaging process using ^* as the next initial sound velocity vector c.

Embodiment 3.
FIG. 11 is a flowchart showing an example of high-precision optical ultrasonic imaging processing according to the third embodiment of the present invention. In the third embodiment, the time is changed by changing the number T of time slots per predetermined unit time (corresponding to the reciprocal of the sampling frequency) from a small number to a large number such as T = 50, T = 200, and T = 1000. The resolution is changed in ascending order, the optical ultrasonic imaging process (S100) according to the first embodiment is executed, and the sound velocity vector c ^* obtained in each process is used as the next initial sound velocity vector c for imaging. It is characterized by performing processing.

In FIG. 11, in step S21, the number of time slots T = 50, the optical ultrasonic imaging process of step S100 is executed in step S22, and the sound velocity vector c ^* is obtained. Next, in step S23, the sound velocity vector c ^* obtained above is used as the initial sound velocity vector for the next imaging process.

Next, in step S24, with the number of time slots T = 200, the optical ultrasonic imaging process of step S100 is executed in step S25 to obtain the sound velocity vector c ^* . Next, in step S16, the sound velocity vector c ^* obtained above is used as the initial sound velocity vector for the next imaging process.

Next, in step S25, the number of time slots T = 1000, the optical ultrasonic imaging process of step S100 is executed in step S26, and the result data (u, c _* ) including the sound velocity vector c ^* is obtained and displayed. Output to 15 for visualization and display.

As described above, according to the third embodiment, the time resolution is changed in ascending order, the optical ultrasonic imaging process (S100) according to the first embodiment is executed, and the sound velocity vector c obtained in each process is performed. By performing the imaging process using ^* as the next initial sound velocity vector c, the accuracy of the imaged image can be improved.

In addition, you may combine Embodiments 2 and 3. That is, at least one of the spatial resolution and the temporal resolution is changed in ascending order, the optical ultrasonic imaging process (S100) according to the first embodiment is executed, and the sound velocity vector c ^* obtained in each process is obtained as follows. Imaging processing may be performed using it as the initial sound velocity vector c.

Embodiment 4.
FIG. 12 is a flowchart showing an image processing using the sound velocity initial value atlas according to the fourth embodiment of the present invention. The imaging process according to the fourth embodiment is characterized in that the imaging accuracy is improved by using the sound velocity initial value atlas. Here, the sound velocity initial value atlas is a two-dimensional or three-dimensional sound velocity atlas that is created and held in advance and is used to determine the initial sound velocity (standard value) in the first embodiment, and the standard sound velocity area is indicated. Refers to the map that has been created.

FIG. 13 is a schematic diagram showing a sound velocity initial value atlas used in the imaging process of FIG. As is clear from FIG. 13, the initial value of sound velocity is different in each region of fat and muscle in the inspection object 1. Here, fat and muscle are examples of regions where the speeds of sound differ from each other.

In step S31 of FIG. 12, first, the sound velocity is fixed at 1500 m / s, and the image data is reconstructed by using a known back projection method (or CS method) (FIG. 14A). Here, if the speed of sound is different, the image quality is affected, but it is considered that the influence is small in knowing the rough location of the tissue. Next, in step S32, the image data obtained in step S31 is subjected to known processing such as smoothing and expansion and known image processing such as threshold processing, and the body containing blood vessels is obtained from the blood vessel distribution. The tissue area is estimated (Fig. 14B). In step S33, in step S32, the body tissue region estimated by using a known alignment method (see, for example, Non-Patent Document 20) and the initial sound velocity atlas are aligned (FIG. 14C). The initial sound velocity map (FIG. 14D) is calculated. Specifically, for example, the body tissue region (inside the alternate long and short dash line) obtained in step S22 is aligned with respect to the atlas region (inside the dotted line) by the following method.
(1) For example, by affine transformation, the body tissue region is aligned with the atlas region by translation, scale, rotation, etc. (see, for example, Non-Patent Document 18).
(2) For example, the body tissue region is aligned with the atlas region by non-rigid body alignment (see, for example, Non-Patent Document 19) with a grid point as a control point.

Next, in step S34, regions having different sound velocities and an initial sound velocity map are generated by using the initial sound velocity map after alignment obtained in step S33. Further, in step S35, the sound velocity map and the reconstructed image data are estimated by executing the optical ultrasonic imaging process (S100) using the initial sound velocity map generated in step S34 as the initial sound velocity.

As described above, according to the present embodiment, by using the sound velocity initial value atlas and aligning the image in combination with the first embodiment, it is possible to image with higher accuracy than the method according to the first embodiment.

Each step of the imaging process or method according to each of the above embodiments can be configured as a computer-executable program. The program can be stored on a computer-readable recording medium such as a CD-ROM or DVD-ROM.

As described in detail above, according to the present invention, it is possible to provide a medical imaging device or the like capable of imaging an object of a living tissue such as a human body with higher accuracy than the prior art.

1 ... Inspection target,
10 ... Optical ultrasonic imaging device,
11 ... Controller,
12 ... Signal receiving circuit,
13 ... Data memory,
14 ... Operation unit,
15 ... Display,
16 ... Image processing unit,
20 ... Optical ultrasonic light source and detection device,
21 ... Sensor support,
22 ... Sensor,
23 ... Mesh board,
24 ... water,
25 ... Laser light source.

Claims (14)

By the optical ultrasonic imaging method, an object is irradiated with light, ultrasonic waves generated from the object are detected as a plurality of detection signals by a plurality of sensors, and from a plurality of detection signals showing a sound pressure spatial distribution. In an optical ultrasonic imaging device provided with a control means for performing an imaging process for imaging the object based on the sound pressure data.
The control means is an objective function relating to the initial sound pressure vector when the multiplication result obtained by multiplying the projection matrix of the sound velocity data by the initial sound pressure vector is used as the sound pressure data based on the predetermined initial sound velocity data. The value of the objective function representing the sparseness of the image data of the object or the sparseness of the data underlying the conversion when a predetermined conversion is performed on the image data is a predetermined threshold value. Imaged image data based on the sound velocity vector is calculated by simultaneously estimating the sound pressure vector and the sound velocity vector by changing the sound velocity vector including the sound velocity in a plurality of regions having different sound velocity as follows. An optical ultrasonic imaging device characterized by this. The optical ultrasonic imaging apparatus according to claim 1, wherein the control means calculates imaged image data based on the sound velocity vector so that the value of the objective function is substantially minimized. .. The objective function is the image data of the object, or the L1 norm, L2 norm, or total variation, or a combination thereof, of the data on the basis of the conversion when a predetermined conversion is performed on the image data. The optical ultrasonic imaging apparatus according to claim 1 or 2. The control means is characterized in that the image data of the object in the process of simultaneously estimating the sound pressure vector and the sound velocity vector, or the image data after the imaging is output to the display unit for visualization. The optical ultrasonic imaging apparatus according to any one of 1 to 3. The control means changes at least one of the spatial resolution and the temporal resolution of the optical ultrasonic imaging method in ascending order, executes the imaging process by the optical ultrasonic imaging method, and performs each imaging process. The optical ultrasonic imaging apparatus according to any one of claims 1 to 4, wherein an imaging process is performed using the sound velocity vector obtained in 1) as the next initial sound velocity vector. The object includes first and second regions having different sound velocities from each other.
The control means
(1) Using a predetermined initial sound velocity vector, the object is imaged using a predetermined imaging method, and image data is calculated.
(2) The second region is estimated by performing predetermined image processing on the calculated image data.
(3) The initial sound velocity map is calculated by aligning the estimated second region with the predetermined sound velocity initial value atlas showing the initial sound velocity value of each region of the object by the predetermined alignment method. And
(4) An optical ultrasonic imaging apparatus according to any one of claims 1 to 5, wherein the imaging process according to any one of claims 1 to 5 is performed using the aligned initial sound velocity map. By the optical ultrasonic imaging method, an object is irradiated with light, ultrasonic waves generated from the object are detected as a plurality of detection signals by a plurality of sensors, and from a plurality of detection signals showing a sound pressure spatial distribution. In the optical ultrasonic imaging method for an optical ultrasonic imaging device provided with a control means for performing an imaging process for imaging the object based on the sound pressure data.
When the control means obtains the sound pressure data as the result of multiplying the projection matrix of the sound pressure data by the initial sound pressure vector based on the predetermined initial sound pressure data, the objective function relating to the initial sound pressure vector. The value of the objective function representing the sparseness of the image data of the object or the sparseness of the base data of the conversion when a predetermined conversion is performed on the image data is a predetermined threshold value. Imaged image data based on the sound pressure vector by simultaneously estimating the sound pressure vector and the sound speed vector by changing the sound pressure vector including the sound speeds in a plurality of regions having different sound speeds or the amount of change thereof as follows. An optical ultrasound imaging method comprising the step of calculating. The optical ultrasound according to claim 7, wherein the control means includes a step of calculating an imaged image data based on the sound velocity vector so that the value of the objective function is substantially minimized. Imaging method. The objective function is the image data of the object, or the L1 norm, L2 norm, or total variation, or a combination thereof, of the data on the basis of the conversion when a predetermined conversion is performed on the image data. 7. The optical ultrasonic imaging method according to claim 7 or 8. The control means includes a step of visualizing the image data of the object in the process of simultaneously estimating the sound pressure vector and the sound velocity vector, or the image data after the imaging by outputting the image data to the display unit. The optical ultrasonic imaging method according to any one of claims 7 to 9. The control means changes at least one of the spatial resolution and the temporal resolution of the optical ultrasonic imaging method in ascending order, executes the imaging process by the optical ultrasonic imaging method, and performs each imaging process. The optical ultrasonic imaging method according to any one of claims 7 to 10, wherein an imaging process is performed using the sound velocity vector obtained in 1) as the next initial sound velocity vector. The object includes first and second regions having different sound velocities from each other.
The control means
(1) Using a predetermined initial sound velocity vector, the object is imaged using a predetermined imaging method, and image data is calculated.
(2) The second region is estimated by performing predetermined image processing on the calculated image data.
(3) The initial sound velocity map is calculated by aligning the estimated second region with the predetermined sound velocity initial value atlas showing the initial sound velocity value of each region of the object by the predetermined alignment method. And
(4) An optical ultrasonic imaging method comprising the step of performing the imaging process according to any one of claims 7 to 11 using the aligned initial sound velocity map. A computer-executable control program for an optical ultrasonic imaging apparatus, comprising each step of the optical ultrasonic imaging method according to any one of claims 7 to 12. A computer-readable recording medium comprising storing the control program of the optical ultrasonic imaging apparatus according to claim 13.