Method for imaging body tissue at high speed and robustly from ultrasonography data using wave form inversion

Abstract

To provide a high-speed and robust method for imaging a body tissue.SOLUTION: A method for imaging a body tissue is composed of steps of: (a) acquiring ultrasonography data; (b) generating one wave motion field using a plurality of CPUs; (c) on a waveform recorded at a reception point from the wave motion field, squaring an error between the waveform and an observation waveform acquired at the step (a) and integrating fc of an equation (5) from all of the CPUs; (d) re-dividing a modified gradient amount and a pseudo-Hessian matrix in the x direction in a region in which the number of CPUs is even; (e) performing re-division in the y direction; (f) using a quasi-Newton's method, calculating a step length required for model update of all wave form inversions; and (g) repeating (b) to (f) to obtain imaging of a viscosity Q value in addition to sound wave propagation velocity, where fc is a cost function and a square sum of residuals between observation data (d_obs) at all reception points at all time and synthetic data (d_syn) calculated using an initial model, and Ns, Nr, and Nt are the number of origination points, the number of reception points, and the number of time samples.SELECTED DRAWING: Figure 19

Description

Technical Field

[0001] The present invention relates to full waveform inversion for imaging the inside of the body with the sound wave propagation speed and the attenuation Q value, which is different from the image (echo image) of the conventional reflecting surface (or scattering point) for the data acquired by ultrasonic medical equipment. The content claimed by the present invention is a method for realizing high-speed calculation of numerical operations that are essential for practical application of full waveform inversion to ultrasonic inspection data.

Background Art

[0002] Full waveform inversion is a visualization technique for underground applying the wave theory developed in the field of seismology (Non-Patent Documents 1 and 2). In principle, this full waveform inversion can be applied to the data acquired by ultrasonic medical inspection equipment, and several past patents have been published in the form combined with the inspection device (Patent Documents 1 to 6).

[0003] Imaging the inside of the body with ultrasonic waves is desirable in that there is no exposure risk compared to X-ray radiography, CT examination, and mammography using X-rays. As a representative of ultrasonic medical equipment, there is an echo examination device. This can provide images in real time. However, in echo examination, it is necessary to estimate the problem area by a skilled doctor or technician performing repeated examinations.

[0004] In mammography, which is a representative of breast cancer examination, it is difficult to distinguish between dense mammary glands and cancerous tissues in the case of dense breasts. In full waveform inversion, since body tissues can be imaged with two physical quantities, the sound wave propagation speed and the attenuation Q value, it is possible to numerically distinguish between mammary glands and cancerous tissues even in dense breasts. Although MRI has slightly lower resolution than X-ray CT, it has no exposure risk and the computational cost is such that imaging results can be obtained on the day of the examination.

[0005] For full waveform inversion, as described in Patent Document 5, in order to rival MRI, it is necessary to provide real-time images, and in order to provide real-time images, it is necessary to reduce the turnaround time to the order of digits. In addition, when full waveform inversion uses data in a limited frequency band (especially a high frequency band that requires a propagation distance of several wavelengths or more with respect to the size of the target), the resulting initial model dependency is high, and it is difficult to ensure robustness. Full waveform inversion in the time domain significantly improves robustness compared to that in the frequency domain, but there are problems with the calculation time and it has not been fully studied. Due to the above problems, full waveform inversion has not yet been put into practical use as an ultrasonic tissue imaging technique at present.

Prior Art Documents

Patent Documents

[0006]

Patent Document 1

Patent Document 2

Patent Document 3

Patent Document 4

Patent Document 5

Patent Document 6

Non-Patent Documents

[0007]

Non-Patent Document 1

Non - Patent Document 2

Non - Patent Document 3

Non - Patent Document 4

Non - Patent Document 5

Non - Patent Document 6

[0008] An object of the present invention is to provide a method for imaging body tissues quickly and robustly by waveform inversion from ultrasonic inspection data in view of the problems of the above prior art and the like. [Means for Solving the Problems]

[0009] First, in order to enhance the robustness of full waveform inversion, a calculation method in the time domain is adopted. In full waveform inversion in the time domain, the image pixel size requires 30 pixels or more per wavelength. This leads to an increase in calculation time but has the advantage of higher image resolution. Instead of the acoustic wave equation cited in the prior patent, a viscous acoustic wave equation is adopted. In addition to the acoustic wave propagation velocity, the viscous Q value can be provided as an imaging result.

[0010] Next, in order to speed up the calculations, parallel calculations are performed on a large number of computers. Since each individual computer may be composed of a general-purpose CPU, the cost is low even when using a large number of computers simultaneously. However, since full waveform inversion requires an enormous amount of memory space, it is necessary to divide the memory. That is, the means for speeding up the calculation of full waveform inversion is summarized as incorporating a method for dividing the memory space into full waveform inversion as follows.

[0011] (1) When solving the wave equation (both forward and inverse propagation), the memory is divided at equal intervals in the x and y directions. In the divided memory space, the wave equation for the initial model is solved using the finite difference method. At each time step of the finite difference method, data exchange is performed at the ends of the divided regions. At the same time, after obtaining the cost function, which is the error from the observed data, in each region, it is collected from the divided memory regions and added together. (2) Calculate the correction amount gradient (Gradient) and the pseudo Hessian matrix as they are in the regions divided in (1). (3) Apply Gaussian smoothing in the y direction to the correction amount gradient and the pseudo Hessian matrix. At this time, divide only at equal intervals in the x direction. (4) Subsequently, apply Gaussian smoothing in the x direction to (3). At this time, divide only at equal intervals in the y direction. (5) While maintaining the memory regions divided at equal intervals in the y direction of (4), correct the correction amount gradient by the pseudo Hessian matrix. (6) While maintaining the division at equal intervals in the y direction of (4), calculate the norm of the correction amount gradient, and collect and add it together into one from the divided memory regions. (7) Calculate the direction vector of the correction amount from the norm obtained in (5), and multiply it by the correction amount gradient to obtain the correction amount for the initial model. Obtain the updated model using the obtained correction amount as the initial model. (8) Return the memory divided regions to the form of (1), and repeat (1) to (7). (9) When the cost function converges below a certain threshold value or the number of repetitions reaches the set upper limit number, end the update operation and use it as the final imaging result.

[0012] That is, the method for imaging body tissue by waveform inversion from ultrasonic inspection data includes: (a) obtaining ultrasonic inspection data; (b) generating one wave field with a plurality of computers (CPUs); (c) integrating, from all CPUs, the squared error between the waveform recorded at the receiving point from the wave field and the observed waveform obtained in step (a), and fc in the following formula (5); (d) re-dividing the correction gradient amount and the pseudo-Hessian matrix in the x-direction into regions with an equal number of CPUs; (e) re-dividing the regions with an equal number of CPUs in the x-direction into regions with an equal number of CPUs in the y-direction; (f) calculating, by the quasi-Newton method, the step length required for model updating of full waveform inversion while keeping the regions divided in the x-direction; and (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the sound wave propagation speed. [Formula 5] TIFF2025090473000002.tif12145Here, fc (scalar quantity) is the cost function, which is the sum of the squared residuals of the observed data (d_obs) and the synthetic data (d_syn) calculated with the initial model at all times for all receiving points. Ns, Nr, and Nt are the number of transmitting points, the number of receiving points, and the number of time samples, respectively. [Advantages of the Invention]

[0013] According to the present invention, a method for imaging body tissue quickly and robustly by waveform inversion from ultrasonic inspection data can be provided. [Brief Description of the Drawings]

[0014]

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15

Figure 16

Figure 17

Figure 18

Figure 19

Figure 20

Figure 21

Figure 22

Figure 23

Mode for Carrying Out the Invention

[0015] For the purpose of explaining the mode, numerical simulations are performed on the numerical model and the arrangement of piezoelectric elements shown in FIG. 1. FIG. 2 shows the physical quantities of body tissues based on experiments. The sound wave propagation speeds of water and tumors are close, but there is a large difference in the viscous Q value. This suggests that when a tumor exists in a place with a lot of moisture in the body tissue (blood vessels or mammary glands), it is difficult to detect the tumor only by imaging the sound wave propagation speed, and the imaging result of the viscous Q value plays a major role in tumor detection.

[0016] The ultrasonic signal shown in FIG. 3 is an ultrasonic signal with a narrow band of about 2 octaves or less that can be generated by a general-purpose piezoelectric element (in the case of FIG. 2, it is in the range of approximately 200 kHz to 800 kHz, but any frequency can be selected for the band).

[0017] In FIG. 1, it is assumed that the geometry surrounds the body tissue with 256 piezoelectric elements as one segment group and 8 segments (2048 piezoelectric elements in total in a clockwise direction). From a practical perspective, there is a certain gap between the segments, and it is assumed that data cannot be acquired in this gap section (radius 11.7 cm, receiver interval 0.336 mm, gap between segments 5.71 mm to 6.38 mm).

[0018] To understand full waveform inversion, consider the case where an ultrasonic signal with the waveform shown in FIG. 5 is transmitted from the piezoelectric element at the position indicated by the arrow (TR1) in FIG. 1.

[0019] The pressure wave transmitted from the piezoelectric element is solved by the finite difference method for the viscous acoustic wave equation expressed by the four simultaneous differential equations of Hooke's law in Equation (1), the memory variable (γ) in Equation (2), and the equations of motion in Equations (3) and (4), to calculate the entire wave field composed of the sound pressure (σ) and the displacement velocities (v x 、v y ). [Equation 1] TIFF2025090473000003.tif11128 [Equation 2] TIFF2025090473000004.tif13149 [Equation 3] TIFF2025090473000005.tif11128 [Equation 4] TIFF2025090473000006.tif11128

[0020] The dot above represents the differentiation in the time direction. ρ and V are the density of the medium and the sound wave propagation velocity, respectively. σ, v x 、v y are the sound pressure, the displacement velocity in the x direction, and the displacement velocity in the y direction, respectively. L is the number of relaxation mechanisms (assumed to be 3). The method of converting from the strict derivation method and the L frequencies (ω l ) defining the viscous Q value to the relaxation time (τ l ) is shown in the inventor's (Non-Patent Document 3).

[0021] The finite difference method is a technique for solving the above differential equations (1) to (4) by directly discretizing (∂x → Δx; ∂t → Δt). To improve the calculation accuracy of the finite difference method, the staggered grid arrangement shown in Figure 4 is adopted. The staggered grid arrangement was devised by Virieux, J., the inventor's supervisor (Non-Patent Document 4). The discretization accuracy adopted is fourth order in the spatial direction and second order in the time direction.

[0022] Figure 5 shows a snapshot of the sound pressure field 64 μs after the transmission time of the forward propagation calculation.

[0023] To speed up the calculations, the calculation area is divided among a plurality of computers at the delimiter positions in Fig. 6 (in the case of Fig. 6, the case of dividing from no division to division by 16 computers is shown, but the number of divisions can be freely set according to the capabilities of the CPUs and the network). Figs. 7 and 8 show diagrams of data exchange between CPUs.

[0024] Fig. 9 shows the time series of the sound pressure received by 2048 piezoelectric elements.

[0025] Next, a method of imaging by full waveform inversion from the simulated observed waveform (Fig. 8) observed by an unknown target (for the sake of explanation, the models in Figs. 1 and 2 are used as the target) will be described. Since an initial model is required for full waveform inversion, an initial model of the sound wave velocity is obtained from the travel time tomography method, which is a general method for obtaining an initial model. Since the viscous Q value cannot be obtained from the travel time tomography method, a model consisting of water and fat excluding the tumor and skin from Fig. 2 is used as the initial model. The density is assumed to be constant at the value of water and will not be solved in the full waveform inversion.

[0026] Fig. 10 shows the initial model obtained by the travel time tomography method. When calculating the "forward propagation wave field in the initial model", synthetic observed data is created at the observation points in Fig. 9 (Fig. 11). Fig. 12 shows the residuals between the observed data (Fig. 9) obtained from the true model and the synthetic observed data (Fig. 11) obtained from the initial model. In the initial model, sound waves in the case where the body tissue does not block between the transmission and reception points can be synthesized in the same way as the observed data. On the other hand, since the initial model has no clear boundaries of body tissues, it can be seen that almost no reflected waves are synthesized.

[0027] The residual is calculated by the sum of the squared residuals of the observed data (d_obs) and the synthetic data (d_syn) calculated with the initial model at all times of all reception points using the following formula 5. [Formula 5] TIFF2025090473000007.tif12145Ns, Nr, and Nt are the number of transmission points, the number of reception points, and the number of time samples, respectively. fc is a scalar quantity and is also called the cost function as it is the objective function to be minimized by the subsequent quasi-Newton method. Since the cost function needs to summarize the values in the divided calculation area, the values of all computers are integrated so that they have the same value on all computers.

[0028] The observed waveform residual is backpropagated (Figure 13 shows a snapshot of the wave field with the residual backpropagated (64 μs after the transmission time)). For the "forward-propagated wave field in the initial model" and the "wave field generated by the backpropagation of the observed waveform residual", the cross-correlation is calculated in the time direction at the imaging points. The zero-lag value of this cross-correlation is the gradient of the model correction amount in full waveform inversion. Instead of storing the forward-propagated wave field in computer memory, the total amount of the correction amount gradient is efficiently obtained by performing the backward calculation of the forward-propagated wave field shown in the inventor's Non-Patent Document 5. At the same time, the pseudo-Hessian matrix is also obtained.

[0029] In the present invention, imaging of the viscous Q value is performed in addition to the acoustic wave propagation velocity. For this purpose, it is necessary to decompose the total amount of the correction amount gradient into the correction amount gradients of the acoustic wave propagation velocity and the viscous Q value. This decomposition method is not shown in Non-Patent Documents 1 and 2. Therefore, the inventor showed the method in Non-Patent Document 6.

[0030] The above-described forward propagation, backward propagation, and calculation of the correction gradient amount are shown together with the following comparative examples to explain the technical feature of performing full waveform inversion in the time domain instead of the frequency domain (which many prior arts use).

[0031] Once the correction amount gradients of the acoustic wave propagation velocity and the viscous Q value are obtained, a Gaussian filter (smoothing) is applied to each correction amount gradient and the pseudo-Hessian matrix. Since Gaussian smoothing cannot be calculated within individual computers while keeping the divided regions intact, the calculation region is divided in the x direction and Gaussian smoothing is performed in the y direction. Next, the calculation region is divided in the y direction and Gaussian smoothing is performed in the x direction (Fig. 15). Due to the characteristics of Gaussian smoothing, the result is the same even if this order (of x and y) is reversed.

[0032] After applying the Gaussian smoothing, the correction amount gradient is corrected by dividing the correction amount gradient (G) by the inverse matrices (H -1 ) of the two pseudo-Hessian matrices for the acoustic wave propagation velocity and the viscous Q value. This correction operation using the pseudo-Hessian matrix can be said to enhance the convergence to the solution of the correction amount gradient without explicitly solving the inverse matrix of the Hessian matrix, which requires a huge amount of computation, by using only the diagonal components of the Hessian matrix that dominate the influence on the correction amount gradient (Non-Patent Document 7). However, it only enables numerical calculation and is theoretically incomplete and empirically works well. Also, an effect of suppressing the generation of locally incorrect images due to geometric inhomogeneity can be expected.

[0033] 〔Equation 6〕 Δm = -H -1 G (6) By adding the correction amount (Δm) of the above Equation 6 thus obtained to the initial model, a corrected model can be obtained. However, since this correction method is a linear approximation, the final imaging is obtained by repeating the calculation from several tens to several hundreds of times with the corrected model. When correcting the model, an operation to accelerate convergence is generally performed by multiplying the correction amount gradient by the step length (α). The calculation for updating from the k-th model (m k ) to the (k + 1)-th model (m k+1 ) is represented by the following Equation 7. 〔Equation 7〕 m k+1 = m k + αΔm k (7) There is also an idea of obtaining the step length from a cost function (evaluated only at the receiving point positions) (adopted by Patent Document 5), but originally it is better to evaluate it for the entire system, so the product (scalar) of the correction amount gradient and the correction amount in the following (8) after correction divide the square (scalar) of the correction amount gradient in the following (9) by TIFF2025090473000008.tif9128 and use the resulting value. i is the grid number and N is the total number of grids. The cost function is used only for determining whether convergence has occurred. So far, it is the steepest descent method, and α guides the objective function to the minimum solution in the tangent direction. However, the correction amount gradient calculated in the time domain has few local errors, and by adopting the quasi-Newton method called the LBFGS (Limited-memory Broyden Fletcher Goldfarb Shanno) method and correcting α in the tangent direction of the quadratic function approximation of the objective function, it is possible to further accelerate convergence.

[0034] However, when the iteration count k = 1, since there is no history of the cost function, the model is updated by the steepest descent method. Furthermore, when the iteration count k = 2, the model is updated by the conjugate gradient method. As a result, for k = 3 and later, at least two cost functions and correction amounts required for the LBFGS method can be obtained.

[0035] The divided area for each computer when performing iterative inversion can remain divided in the y direction after the Gaussian filter. In the calculation of model update applying this inversion technique, the sum of the squares of the correction gradient amounts and the sum of the products of the correction gradient amounts and the corrected correction amounts are compared with their histories for execution. Therefore, it is necessary to integrate the L2-norm from the areas divided in the y direction so that all computers have the same value.

[0036] The imaging result obtained as described above is shown in FIG. 16.

[0037] Plot the cost function against the number of iterations in Figure 17. It can be seen that the cost function decreases as the number of iterations increases.

[0038] In Figure 18, arrange the updated models at the 10th, 20th, and 30th iterations. The focus of the tumor can be captured within 10 iterations, but it takes 20 - 30 iterations for the skin to come into focus.

[0039] The acoustic wave propagation speed is imaged accurately. There is a certain error in the viscous Q value. Part of the cause of this error is that the density model is kept constant. However, most of the problem is that, as can be seen by looking around the skin, the image is vibrating, which means that the acoustic wavelength is too short relative to the size of the target. Nevertheless, the relative magnitude relationship of the viscous Q value image preserves the true model and is accurate enough for diagnosis as an image.

[0040] Summarize the flowchart in Figure 19. The method for decomposing the calculation region of the present invention is illustrated along the flowchart of full waveform inversion.

[0041] Fig. 20 shows the results of benchmark tests on calculation time and memory. In the benchmark test, the number of pixels in the x direction is 3005, the number of pixels in the y direction is 3005, the number of samples in the time direction is 5625, and the number of transmission points is 24. The computer consists of 2 Intel XEON (E5-2690 v4 with a clock frequency of 2.6 GHz, 14 cores) CPUs and 256 GB of memory to form one server, and a PC cluster with 24 servers connected by LAN cables is used. For parallel numbers of 24, 48, 96, 192, and 384, the values normalized with respect to the parallel number of 24 are displayed. They respectively mean that the number of domain partitions during finite difference method calculations is 1, 2, 4, 8, and 16. When the number of domain partitions is 2, the calculation time is reduced to less than half, and the improvement in the efficiency within the computer due to the reduction of memory exceeds the loss due to communication. When the number of domain partitions is 2 or 4, the length of the sides of the calculation area that requires communication is homogeneous in all computers, but for 8 and 16, the length of the sides for communication is longer for the computers calculating the central part. It can be seen that the loss due to this inhomogeneity of communication volume increases. However, both the calculation time and memory can be reduced approximately proportionally to the parallel number.

[0042] 〔Comparative Example〕 〔Full Waveform Inversion in Frequency Domain〕 The differences between performing full waveform inversion in the time domain and the frequency domain are explained. Using the five frequency samples (251770, 375748, 501633, 625610, 751495 Hz) shown in Fig. 2(b), the calculation of the modified gradient amount was performed for full waveform inversion in the frequency domain.

[0043] Fig. 21 shows the wave field in the frequency (cosine transform) domain at a frequency of 501633 Hz. Since the wave field in the frequency domain is a single harmonic response, it has a striped pattern. As a point to note, strong responses appear outside the body tissue and outside the circle of the transmitter-receiver point arrangement. However, it was decided to perform backpropagation for the recordings of the receiver points within an aperture width of 90° from the transmitter point. This is because when performing backpropagation for the recordings of all receiver points with an aperture width of 180°, the results of the modified gradient amount deteriorated.

[0044] Figure 22 shows the modified gradient amount. In the frequency domain, the modified gradient amount is the real part of the complex product of the forward-propagating wave field and the backward-propagating wave field. However, outside the body tissue, large values also remain outside the circle of the transmitting and receiving point arrangements. This presents a completely different aspect from the modified gradient amount obtained in the time domain of Figure 14. This is the cycle skipping problem that has long been regarded as an issue in full waveform inversion. To reduce this problem, various smoothing and regularization calculations have been devised, but no solution method has been found as far as the inventor has tried.

[0045] The modified gradient amount also increases at a certain period within the body tissue, showing a concentric circle-like pattern. This means the period at which the phases of the five frequency samples used are aligned. Unfortunately, even if the frequency samples are changed within the above range, due to the very short wavelength (about 3 mm) relative to the target size (16 cm), a period at which the phases are aligned (reinforced or weakened) will necessarily occur. To avoid this, all frequency samples can be calculated, but trying this is the same as full waveform inversion in the time domain and is extremely heavy compared to the calculations in the time domain. The cycle skipping problem can be solved by lowering the frequency band of the signal, but for that, it is necessary to generate sound waves that penetrate at 1 to 10 wavelengths relative to the target size. To mechanically generate low frequencies, it is necessary to increase the size of the piezoelectric element, which goes against the significance of low-cost inspection equipment, which is an advantage of ultrasonic medical devices.

[0046] Also, there is a secondary wave field generated by the gap between segments. It can also be seen that the modified gradient amount has a large reaction to this artificially generated noise. Therefore, although patent documents state that implementation is possible in either the frequency domain or the time domain, with regard to the frequency domain, it is insufficient in that it is only possible if calculations can be performed for all frequencies. For example, in Patent Document 4, full waveform inversion is referred to as characterization of the boundary region and is only used to enhance X-ray CT and travel time tomography images (initial models), and those that cannot be imaged to a certain extent with the initial model will not be imaged by full waveform inversion. This implicitly indicates that the robustness of full waveform inversion in the frequency domain is low. Therefore, the applicant claims that it is only possible to robustly obtain imaging results in the time domain (Fig. 23). Furthermore, it is considered that by obtaining an image of the viscous Q value in addition to the sound wave propagation velocity, the probability of detecting a tumor can be improved.

Industrial Applicability

[0047] According to the present invention, since it is possible to provide a method for imaging body tissues quickly and robustly by waveform inversion from ultrasonic inspection data, it is possible to image body tissues safer, at lower cost, and with higher resolution than X-ray CT or MRI of conventional inspection systems. In particular, conventionally, by effectively performing parallel calculations on a large number of CPUs for the problem of the calculation time of waveform inversion in the time domain, it can be said that there is a great possibility of industrial use in that analysis results can be obtained in a calculation time comparable to that of conventional inspection systems.

Explanation of Signs

[0048] T / R1 Transmission point of snapshot (Fig. 5), first reception point T / R Transmission and reception point (A) Tumor (B) Fat (C) Skin (D) Water V Sound wave propagation velocity Q Viscous value ρ Density of medium

Claims

1. An imaging method of body tissues by waveform inversion from ultrasonic inspection data, comprising: (a) obtaining ultrasonic inspection data; (b) generating one wave field with a plurality of computers (CPUs); (c) for the waveform recorded at the receiving point from the wave field, integrating the squared error between the observed waveform obtained in step (a) and fc in the following formula (5) from all CPUs; (d) re-dividing the correction gradient amount and the pseudo-Hessian matrix in the x direction into regions with an equal number of CPUs; (e) re-dividing the regions with an equal number of CPUs in the x direction into regions with an equal number of CPUs in the y direction; (f) calculating, by the quasi-Newton method, the step length required for model updating of full waveform inversion while keeping the regions divided in the x direction; and (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the sound wave propagation velocity, a fast and robust imaging method of body tissues. [Formula 5] Here, fc (scalar quantity) is a cost function, which is the sum of the squared residuals of the observed data (d_obs) and the synthetic data (d_syn) calculated with the initial model at all times of all receiving points. Ns, Nr, and Nt are the number of transmitting points, the number of receiving points, and the number of time samples, respectively.

2. In step (b), data exchange is performed at each time step, the number of divisions of the divided region is set according to the computer environment (number of CPUs, memory), the divided region is applied to all three wave field calculations of the forward propagation, backward propagation, and backward calculation of the forward propagation of full waveform inversion, and the correction amount gradient and the pseudo-Hessian matrix are obtained while keeping the divided region, a fast and robust imaging method of body tissues according to claim 1.

3. After the step (d), a Gaussian filter is applied in the y direction, and after further re-dividing from an equal area in the x direction into an equal area of the number of CPUs in the y direction, a Gaussian filter is applied in the y direction. The method for fast and robust imaging of body tissues according to claim 1, characterized in that.

4. The sum of squares shown in the following (9) of the correction gradient amount calculated in the step (f), the correction gradient amount, and the sum of products shown in the following (8) of the pseudo-Hessian matrix are integrated from all CPUs. The method for fast and robust imaging of body tissues according to claim 3, characterized in that. 〔8〕 〔9〕 Here, i is the grid number and N is the total number of grids, H －１ is the inverse matrix of the pseudo-Hessian matrix and G is the correction amount gradient.

5. The ultrasonic inspection data is ultrasonic breast cancer inspection data, The body tissue includes skin, breast, and cancerous tissue. The method for fast and robust imaging of body tissues according to claim 4, characterized in that.

Images