Ultrasonic Tomographic Technique and Its Applications

Abstract

X-ray tomography and magnetic resonance imaging (MRI) are excellent techniques for non-destructive or non-invasive inspections, however, they have shotcomings including the expensive cost in both the devices themselves and their protection facilities, the harmful side effects of the X-rays to human bodies and to the environment. In view of this argument, it is necessary to develop new, inexpensive, safe and reliable tomographic techniques, especially in medical imaging and non-destructive inspections. There are new tomographic techniques under development such as optical tomography, photo-acoustic tomography, ultrasonic tomography and so on, from which we take ultrasonic tomography as the topic in this paper. We introduce a review of the known ultrasonic tomographic techniques and discuss their future development.

1. Introduction

In this paper, we review the known techniques in the development of ultrasonic tomography and discuss their future development. As for reliable and minute tomographic techniques, the X-ray computerized tomography (X-ray CT) and the magnetic resonance imaging (MRI) are well known and frequently applied, especially for non-invasive inspections in medical imaging. They are superior techniques in non-invasive and non-destructive inspections, however, they have deficiencies such as the expense of both the devices themselves and their protective facilities, the harmful side effects of the X-rays to human bodies, the unsuitability of an MRI for a non-destructive inspection of those containing steel inside their bodies and so on. In order to avoid such shortcomings, tomographic techniques other than X-ray CT and MRI are under investigation; for example, electrical impedance tomography, acoustic tomography, optical tomography, photo-acoustic tomography, ultrasonic tomography and so on.

In the next section, we briefly review tomographic techniques without X-ray nor MRI, among which includes ultrasonic tomography which we shall be interested in.

In the third section, we first study why G. N. Hounsfield was successful in practicalizing CT. We briefly review his idea and give a theoretical proof why his idea worked well, based on which, we propose a new algorithms for ultrasonic CT for its further development.

In the last section, we summarize the conclusion in this paper.

2. New Tomographic Techniques

As mentioned in the previous section, new tomographic techniques distinct from X-ray and MRI are under investigation from the viewpoint of both theory and practice. In this section, we review such new tomographic techniques, among which, we shall be interested in ultrasonic tomography (or ultrasonic computerized tomography (USCT)). We also introduce why we shall be interested in ultrasonic tomography at the end of this section.

2.1. Acoustic Tomography

The history of acoustic tomography is old. It is based on inverse problems of the wave equation as Equation (1);

$$({\partial }_t^2-C\left(x\right)\Delta )u(x,t)=0,$$

(1)

where $C\left(x\right)$ is the velocity of the wave at the point $x\in \mathsf{\Omega }$ with some domain $\mathsf{\Omega }\subset {\mathbb{R}}^3$, and $t>0$. In this problem, by controlling and observing each of the two boundary data as Equation (2)

$$u{|}_{\partial \mathsf{\Omega }}$$

(2)

or as Equation (3)

$$\gamma \frac{\partial u}{\partial \nu }{|}_{\partial \mathsf{\Omega }},$$

(3)

we try to reconstruct the function $C\left(x\right)$. Mathematical research in acoustic tomography became the origin of the study of what are called “Neumann-to-Dirichlet” and “Dirichlet-to-Neumann” maps. Acoustic tomography is an inverse problem to reconstruct the wave velocity at each space point from the control and observation of the boundary data. Consult the second section in [1] for the review of the acoustic tomography. This idea is applied for geophysical analysis [2], ocean acoustic tomography [3] and others.

2.2. Electrical Impedance Tomography (EIT)

In EIT, we study inverse problems of the equation as Equation (4)

$$\nabla ·\gamma (x,\omega )\nabla u\left(x\right)=0,$$

(4)

where $u\left(x\right)$ is the electric potential at the point $x\in \mathsf{\Omega }$ in the body $\mathsf{\Omega }\subset {\mathbb{R}}^3$, $\gamma (x,\omega )=\sigma (x,\omega )+i\omega ϵ(x,\omega )$ is the admittivity, $\sigma (x,\omega )$ is the electric conductivity, $ϵ(x,\omega )$ is the electric permittivity and $\omega$ is the angular frequency of the applied current. In Equation (4) if we control the input current density Equation (3) where $\partial \mathsf{\Omega }$ is the boundary of $\mathsf{\Omega }$ and $\nu$ is the outward normal vector on $\partial \mathsf{\Omega }$, and we observe the voltage Equation (2) on the boundary then we obtain a Neumann-to-Dirichlet type problem. Conversely, if we control Equation (2) and observe Equation (3) then we obtain a Dirichlet-to-Neumann type problem. EIT is an inverse problem to reconstruct the admittivity at each space point and a frequency from Neumann-to-Dirichlet or Dirichlet-to-Neumann type data. Consult [4] for the review of EIT technique and [5], for example, for an application of EIT for medical imaging.

2.3. Optical Tomography

In optical tomography, we project the infrared lazer beams to an object and observe their advection through and scattering by the object. From the observation on the boundary of the object, we try to reconstruct the interior information, which is formulated an an inverse problems in the transport equation as Equation (5)

$${\partial }_{t}u(x,\theta ,t)+v·\theta \nabla u(x,\theta ,t)-{\mu }_{a}u(x,\theta ,t)={\mu }_s{\int }_{S^2}\varphi ({\theta }^{\prime },\theta )u(x,{\theta }^{\prime },t)d{\theta }^{\prime },$$

(5)

where $x,v\in {\mathbb{R}}^3,t\le 0,\theta \in S^2$, v is the velocity vector, $u(x,\theta ,t)$ is the density of photon at $(x,t)$ in the direction $\theta \in S^2$, ${\mu }_a$ is the absorption constant and ${\mu }_s$ is the scattering constant. The right hand side of Equation (5) represents the effect of scattering, the last term in the right hand side represents the effect of absorption and the second term in the right hand side represents the advection effect. In general it being difficult to exactly know the scattering kernel $\varphi ({\theta }^{\prime },\theta )$, suitable a priori informations in accordance with the phenomena are known. In the study of optical tomography, reconstruction of any of $v,{\mu }_a,{\mu }_s$, with some a priori assumption on the others, can be of interest. For the introductive and general theory of the optical tomography, consult [6]. If the advection and the absorption effects are small, then it is possible to approximate the Equation (5) by the diffusion equation as Equation (6)

$${\partial }_{t}u(x,\theta ,t)-\nabla ·D\left(x\right)\nabla u(x,\theta ,t)=0$$

(6)

and reconstruction of the diffusion coefficient $D\left(x\right)$ gives a rough sketch of the interior information, for which consult [7].

2.4. Photoacoustic Tomography (PAT)

PAT is based on the photoacoustic effect. If non-ionized lazer beams are projected to a human body, they go through the skin and reach the interior organs and cells, where lazer energy is absorbed and changed into heat, which projects ultrasound waves. We can observe this ultrasound and change it into the image of the human body. In this imaging process, we solve an inverse problem for the photoacoustic wave equation as Equation (7)

$$({\partial }_t^2-C\left(x\right)\Delta )u(x,t)=\gamma E(x,t),$$

(7)

where $C\left(x\right)$ is the velocity of the wave at the point $x\in \mathsf{\Omega }$ of the human body $\mathsf{\Omega }\subset {\mathbb{R}}^3$, $t\ge 0$, $E(x,t)$ is the amount of absorbed light energy at $(x,t)$ and $\gamma$ is what is called Grüneisen parameter which represents the rate where $E(x,t)$ is changed into ultrasound. It looks that there are a number of unsolved, or even unestablished, mathematical problems for further development of PAT technology. Consult [8] for the survey of PAT.

2.5. Ultrasonic Computerized Tomography (USCT)

Ultrasonic tomography is a kind of an acoustic tomography. In this technique, we apply ultrasound whose frequency differs in accordance with the object for inspection. If the object is concrete structures, we apply ultrasound of the frequency about 50 KHz, while ultrasound of the frequency about 2 MHz is applied for the human body. Though we are able to study the same inverse problems in acoustic tomography, in this paper, let us study other problems than acoustic tomography, which are very special to ultrasonic tomography. For the propagation of the ultrasound, the following property is known.

Property 1.

Let $\mathsf{\Omega }\subset {\mathbb{R}}^3$ be a domain where the object locates, and $C\left(x\right)$ be the propagation speed of the ultrasound at the point $x\in \mathsf{\Omega }$. For $\alpha ,\beta \in \partial \mathsf{\Omega }$, we denote by ${\gamma }_{\alpha ,\beta }$ a route from α to β contained in Ω. The primary wave of the ultrasound which travels from the point $\alpha \in \partial \mathsf{\Omega }$ to the point $\beta \in \partial \mathsf{\Omega }$ takes the route where the travel time is given by Equation (8)

$$\underset{{\gamma }_{\alpha ,\beta }}{min}{\int }_{{\gamma }_{\alpha ,\beta }}1/C\left(x\right)d\gamma ,$$

(8)

This route is called ‘the fastest route’.

For this property consult [9]. By virtue of this property, the mathematical problem for ultrasonic tomography is given as

Problem 1

(Problem to develop USCT).  Let $\mathsf{\Omega }\subset {\mathbb{R}}^3$ be a domain and $C\left(x\right)$ be the propagation speed of the ultrasonic wave at the point $x\in \mathsf{\Omega }$. For $\alpha ,\beta \in \partial \mathsf{\Omega }$, we denote by ${\gamma }_{\alpha ,\beta }$ a route from α to β contained in Ω. In this case, reconstruct $C\left(x\right)(x\in \mathsf{\Omega })$ out of the data as Equation (9)

$$\underset{{\gamma }_{\alpha ,\beta }}{min}{\int }_{{\gamma }_{\alpha ,\beta }}1/C\left(x\right)d\gamma ,$$

(9)

for $\forall \alpha ,\beta \in \partial \mathsf{\Omega }$.

This problem is also posed in [9]. With the mathematical solution to Problem 1 being very difficult and open, we can establish and solve practical problems related to Problem 1. In the paper [9], the authors posed and solved how to non-destructively inspect the concrete cover for reinforcement in reinforced concrete structures.

If the object is small or close to homogeneous, then we can apply the same algorithm with the X-ray CT, that is, we assume that the ultrasound rectilinearly propergates and apply the same algorithm with the X-ray CT to obtain a rough sketch of the distribution of the reciprocal $1/C\left(x\right)$ of the velocity. This idea is successful in medical imaging where the object is very small (cf. [10]). In the next section, we shall extend the this idea for further development of ultrasonic tomography.

2.6. Summary

In this subsection, as the summary of this section, let us compare the tomographic techniques introduced in this section. Before the summary, let us note that there is another tomographic technique called magnetic resonance electrical impedance tomography (MREIT, [11]). This technique is a combination of MRI and EIT, which is not a tomographic technique without X-ray nor MRI, hence omitted to be introduced in this section, for the reasons mentioned in the introduction of this paper.

In acoustic tomography, the device being not so expensive, a number of observation data are required for the reconstruction of the velocity distribution $C\left(x\right)$. We cannot say that the reconstruction of the acoustic tomography is very sharp, for the time being. If a rough sketch of the velocity distribution $C\left(x\right)$ is sufficient or if the object is very big, like an ocean tomography, then this method can play an important role (cf. [3]).

In EIT, mathematical treatment of the problem being same as acoustic tomography, the device is more expensive and the reconstruction is sharper, but more development is still required to function as a sharp medical imaging technique like CT and MRI.

Optical tomography and PAT being new and interesting non-invasive inspection techniques, they should also be more improved as sharp medical imaging techniques.

In many cases, these new techniques are applied to obtain some rough sketches of the object for examination, not for an inspection. Note that by ‘examination’, especially in medical application, we mean the cheaper-running and rougher-sketch-obtaining inspection in order to determine whether a more detailed inspection is necessary or not.

We would like a new technique which can be an alternative to X-ray CT or MRI for non-destructive or non-invasive inspection of the object, not for rough examination. In this view, ultrasonic tomography seems to be possible to be the sharpest imaging technique by virtue of Property 1 and Problem 1, where we can observe as sharp data as the X-ray tomography except that the orbit of the ultrasound is not a line. However, ultrasonic CT cannot be applied for the non-invasive inspection for the human brain, because of Property 1. The ultrasound travels along the skull, where the ultrasonic speed is very fast, and would not reach the interior brain, in view of which, optical tomography is expected to be an alternative to brain CT or brain MRI. While PAT, by its physical nature, is expected to be an alternative of positron emission tomography (PET), where the X-rays are discharged at the expected interior point in the human body. We note that the required device is the cheapest in ultrasonic CT, among these three new tomographic techniques.

In spite of the above demerit, the ultrasonic tomography is possible to be the sharpest imaging technique in the new tomographic techniques introduced in this section, by virtue of Property 1 and Problem 1. This technique can be an alternative of CT and MRI except for the non-invasive inspection of the brain. One of the best merits of the ultrasonic tomography is that the device and its running cost is much cheaper than CT and MRI and that no protective facilities are required. It may enable small clinics and hospitals to have their own CT devices, which can contribute to the development of clinical medicine all over the world. We also note that the cost of the ultrasonic tomography can be the cheapest in the new tomographic techniques introduced in this section. By the inexpensiveness of its cost and by the ease of its application, the ultrasonic tomography can be an effective technique to non-destructively inspect the concrete structures, where it is not suitable to spend too much cost since the concrete itself is very cheap.

In view of the above argument, we shall focus on the ultrasonic tomography and study it in the next section.

3. New Methodology for Ultrasonic CT

It is our main purpose in this paper to propose a new methodology for ultrasonic CT, which shall be discussed in this section.

As mentioned in the previous section, if the object is small or close to a homogeneous one, then, assuming that the ultrasound rectilinearly propagates, we can apply the same algorithm with the X-ray CT. For the time being, the most widely applied CT algorithm is what is called ‘the filtered back projection technique’.

Before discussing a new methodology for USCT, we shortly shortly review the X-ray CT, which shall be the basic idea for our development of a new methodology for USCT. Consider a section of the human body by a plane, where we define the coordinate in order that this plane is given by $\{(x,y,z)\in {\mathbb{R}}^3|z=0\}$. Let $f(x,y)=f(x,y,0)$ be the density of the human body in the plane $\{(x,y,z)\in {\mathbb{R}}^3|z=0\}$. The mathematical problem of CT is introduced as follows.

Problem 2.

Reconstruct the function $f(x,y)$ defined on ${\mathbb{R}}^2$ by its line integrals ${\int }_{l}f(x,y)dl$ along all lines l’s in ${\mathbb{R}}^2$.

Roughly speaking, the problem of X-ray CT is to reconstruct the density distribution $f(x,y)$ for $(x,y)\in {\mathbb{R}}^2$, from its all line integrals ${\int }_{l}f(x,y)dl$.

For the properties of the X-rays, introduction of the mathematical model for CT and the introduction of Problem 2, consult [12,13,14] for example. In Problem 2, it is sufficient to assume $f\in L^1\cap L^2\left({\mathbb{R}}^2\right)$ which is satisfied in the practical CT and to take almost all lines l’s in ${\mathbb{R}}^2$. There happens a phenomenon that the X-ray is attenuated, whose cause is the density $f(x,y)$ of the human body. Therefore, Problem 2 is an inverse problem to reconstruct the cause $f(x,y)$ of the phenomena by the knowledge of its results ${\int }_{l}f(x,y)dl$’s, by which we mean that this problem is an inverse problem to solve the inverse cause-result relation.

3.1. Filtered Back Projection

It is known that the mathematical model of the X-ray CT is given by the Radon transform in the two dimensional Euclidean plane, for which consult [13]. In the filtered back projection, we apply the low-pass filter between the integral and the operator $I^{-1}$ in the following reconstruction formula for the Radon transform, as Equation (10)

$$f\left(x\right)=\frac{1}{4\pi }{I}^{-1}\left({\int }_{S^1}Rf(\theta ,x·\theta )d\theta \right),$$

(10)

where for $x\in {\mathbb{R}}^2$, $\theta \in S^1$ with $S^1$ being the unit circle in ${\mathbb{R}}^2$, $s\in \mathbb{R}$, shown as Equation (11) below

$$Rf(\theta ,s)={\int }_{x·\theta =s}f\left(x\right)dx$$

(11)

is the Radon transform and Equation (12)

$$I^{-1}\phi \left(x\right):=F^{-1}\left(\right|\xi \left|F\phi \left(\xi \right)\right)\left(x\right),$$

(12)

with Equation (13)

$$Ff\left(\xi \right):={\int }_{{\mathbb{R}}^n}{e}^{-ix·\xi }f\left(x\right)dx,F^{-1}\phi \left(x\right):=\frac{1}{{\left(2\pi \right)}^n}{\int }_{{\mathbb{R}}^n}{e}^{ix·\xi }\phi \left(\xi \right)d\xi .$$

(13)

Therefore, the filered back projection given in Equation (10) is Equation (14)

$$f\left(x\right)=\frac{1}{4\pi }{I}^{-1}\mathcal{F}\left({\int }_{S^1}Rf(\theta ,x·\theta )d\theta \right),$$

(14)

where $\mathcal{F}$ is a low-pass filter. It is because that the operator $I^{-1}$ sensitively responds to the errors, some regularization is required, which is realized by a low-pass filter $\mathcal{F}$. Today this technology is applied to almost all medical X-ray CT devices with the low-pass filter $\mathcal{F}$ developed by each CT devise-making company.

In cases where the object is small or close to a homogeneous one, by assuming that the ultrasound rectilinearly propagates, the filtered back projection technique for the X-ray CT is applied for the ultrasonic CT, however, it is only to obtain some rough image by the ultrasonic CT. With this method, further development cannot be expected since the essential orbit of the ultrasound is different from the one of the X-ray.

3.2. The Idea by G.N. Hounsfield

The author claims that, for further development of ultrasonic CT, it is of interest to modify the original idea by G.N. Hounsfield who first practicalized the X-ray CT by developing an device for X-ray CT inspection. In this section, we review his idea, which shall be modified for ultrasonic CT in the Section 3.4. There being few medical CT devices applying G.N. Hounsfield’s idea today, if the object is simple, for example, CT scan for the teeth and so on, then G.N. Hounsfield’s algorithm still works well.

In numerical analysis, we usually directly discretize a mathematical formula obtained by mathematical analysis in order to implement it for practical application, however, the essence of G.N. Hounsfield’s idea is to discretize the mathematical model itself, not a mathematical formula. It is very interesting and worked very well at the initial stage of the practicalization of CT. We can assume that $\mathrm{supp}f$ is compact and $f\in L^1\left({\mathbb{R}}^2\right)$ since $f(x,y)$, $(x,y)\in {\mathbb{R}}^2$ is a two-dimensional section of the density distribution of the human body. Cover $\mathrm{supp}f$ with n squares $c_1,c_2,\cdots ,c_n$, whose sides are of the same length and parallel to the x- or y-axis and any pair of whose interiors are mutually disjoint. We approximate the function f by a function $g(x,y)$ defined as Equation (15)

$$g(x,y)=\sum _{j=1}^n{x}_j{\chi }_{c_j}(x,y),$$

(15)

where ${\chi }_{c_j}(x,y)$ is the characteristic function of the square $c_j$ and $x_j$ are unknowns. For example, if we suppose the unknown $x_j$ as the integral mean of $f(x,y)$ in the square $c_j$ as Equation (16)

$$x_j=\frac{1}{|c_j|}{\int }_{c_j}f(x,y)dxdy,$$

(16)

then it seems easy to understand G.N. Hounsfield’s idea, where $|c_j|$ is the area of the square $c_j$. The function $g(x,y)$ is called a pixel function.

For example, in Figure 1, the original image (left) of the section of a human body is divided into the sum of squares (right) and in each square the density is assumed to be a constant.

G.N. Hounsfield tried to reconstruct an approximation of the pixel function $g(x,y)$, neither $g(x,y)$ itself nor $f(x,y)$ itself, from the observed data necessarily containing errors in various senses. Since the best we can hope is to obtain an approximation of f, not to reconstruct f itself, it is very flexible to approximate f by a suitably simple function g so that the problem would be simplified. Assume that m X-rays, $I_1,I_2,\cdots ,I_m$, are projected to the human body. We assume that for $i=1,2,\cdots ,m$, the strength of $I_i$ before the projection is $I_i^0$ and after the projection is $I_i^1$. By $l_1,l_2,\cdots ,l_m$, we denote the lines where the X-ray $I_1,I_2,\cdots ,I_m$ rectilinearly propagates, respectively and by $a_{ij},i=1,2,\cdots ,m,j=1,2,\cdots ,n$, we denote the length of $l_i\cap c_j$. Letting $s_i:=log{I}_i^0-log{l}_i^1,i=1,2,\cdots ,m$ together with Equation (15), the problem we have to solve turns out to be the following one.

Problem 3.

Solve the following system of linear equations in $x_1,x_2,\cdots ,x_n$ as Equation (17).

$$\left\{\begin{array}{cc}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n& =s_1,\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n& =s_2,\\ \cdots \cdots \cdots \\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n& =s_m,\end{array}\right.$$

(17)

or equivalently as Equation (18)

$$A\mathit{x}=\mathit{s}.$$

(18)

In the medical CT, we have to treat a huge system like $m\ge 9\times {10}^6,n\ge 3\times {10}^4$. The fact $m>>n$ and the effect of the errors we mentioned above, the overdetermined system Equation (17) must have no solution. Although the system Equation (17) has no solution, the density distribution $f(x,y)$ of the human body exists with no doubt, which must be approximately solved. Therefore, we have to solve the following problem.

Problem 4.

Solve the overdetermined system of linear equations in Equation (17) with no solution.

This is the problem G.N. Hounsfield tried to solve. His idea is what is now called algebraic reconstruction technique (ART) or Kaczmark’s method. He discovered this idea independently without knowing ART nor Kaczmarz’s method (for these methods, consult [13]). In the rest of this section, we shall introduce the G.N. Hounsfield’s solution of Problem 4.

We first note that each line as Equation (19)

$$a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n=s_i$$

(19)

in Equation (17) is a equation of the hyperplane $H_i\subset {\mathbb{R}}^n$ whose normal vector is ${\mathit{a}}_i:=(a_{i1},a_{i2},\cdots ,a_{in})$. Take any ${\mathit{x}}_0\in {\mathbb{R}}^n$ and define the sequence $\left\{{\mathit{x}}_k\right\}$ by Equation (20)

$${\mathit{x}}_k:=P{\mathit{x}}_{k-1}\equiv P_m{P}_{m-1}\cdots P_1{\mathit{x}}_{k-1},k=1,2,\cdots$$

(20)

where Equation (21)

$$P_i\mathit{x}:=\mathit{x}-\frac{{\mathit{a}}_i·\mathit{x}-s_i}{\parallel {\mathit{a}}_i\parallel }\frac{{\mathit{a}}_i}{\parallel {\mathit{a}}_i\parallel }$$

(21)

is the orthogonal projection of $\mathit{x}$ onto the hyperplane $H_i$. If the system Equation (17) has the unique solution, then it is easily proved that Equation (22)

$$\underset{k\to \infty }{lim}{\mathit{x}}_k$$

(22)

converges to the unique solution of Equation (17), which is called ART or Kaczmark’s method (confer Figure 2 for its image).

Of course, the system has no solution in practice, because of errors in approximation of the original function f by the pixel function g, errors in observation, numerical rounding errors by computers, some noises and so on, all of which are small. Therefore, it can be said that the system Equation (17) is “close to have the unique solution”. Therefore G.N. Hounsfield considered that Equation (22) would give an approximate solution of Equation (17), which worked very well in practice. In fact, we find that the iteration method by G.N. Hounsfield almost converges in a few steps, $k=3$ or 4 in Equation (22), by numerical experiments. There being a number of works relating the justification why his idea works well, consult [15,16] for example, it seems that there in no direct proof to justify his idea. In the next subsection, we shall give an approach to give a theoretical proof why G.N. Hounsfield’s idea worked well.

3.3. A Theoretical Justification of G.N. Hounsfield’s Idea

We give a justification why G.N. Hounsfield’s idea worked well in practice. His iteration method is justified in terms of the least square solutions. If a linear system Equation (17) has a solution $\mathit{x}$, then it holds $A\mathit{x}-\mathit{s}=\mathbf{0}$, in view of which the following the idea of the least square solution is defined.

Definition 1.

A vector$\mathit{x}\in {\mathbb{R}}^n$is called a least square solution of Equation (17) (or Equation (18)) if and only if it minimizes the norm as Equation (23)

$$\parallel A\mathit{x}-\mathit{s}\parallel ,$$

(23)

where$\parallel \mathit{y}\parallel :=\sqrt{y_1^2+\cdots +y_m^2}$for$\mathit{y}\in {\mathbb{R}}^m$.

We note that there are many researches to give least square solutions to the overdetermined system Equation (17), consult [15,16] for example, our theory is different from them in the following points.

• We a priori know that the least square solution to the overdetermined system Equation (17) is unique (Theorem 1 below).
• We shall discuss not only the solution to Equation (17) but also the limit where the size of the all pixels tends to 0.
• We know by experience that, even for fixed number of pixels (even if the number n is fixed), the more observation X-ray data we have (the larger the number m is), the better resolution we can obtain, which shall be discussed in a stochastic way (Theorem 2 below).

As for the least square solutions, there hold the following propositions.

Proposition 1.

For $A\in M_{mn}\left(\mathbb{R}\right)$, the following conditions are equivalent.

(i)

There exists at least one least square solution of Equation (17) for any $\mathit{s}\in {\mathbb{R}}^m$.

(ii)

Furthermore, if rank $\left(A\right)=n$ then the solution of Equation (17) is unique.

Proposition 2.

For $A\in M_{mn}\left(\mathbb{R}\right)$, the following conditions are equivalent.

(i)

$\mathit{x}\in {\mathbb{R}}^n$ is a least square solution of Equation (17).

(ii)

There holds the following Equation (24).

$${}^{t}AA\mathit{x}-{}^{t}A\mathit{s}=\mathbf{0}.$$

(24)

The proofs of Propositions 1 and 2 being exercises of linear algebras, we would not mention them. In the practical CT, the rank of the coefficient matrix A in Equation (17) must be n. Therefore, we have the following theorem.

Theorem 1.

In the practical CT, there always exists the unique least square solution to the system Equation (17) of the linear equations.

In the following, we assume the following conditions.

(I)

$f(x,y)\in L^1\cap L^2\left({\mathbb{R}}^2\right)$.

(II)

For any line $l\subset {\mathbb{R}}^2$, ${f|}_l\in L^1\left(l\right)$ and is piecewisely continuous.

Note that these conditions, (I) and (II), are satisfied by the density distribution of the human body. We also assume that

(III)

$m>>n$.

With these assumptions, (I), (II) and (III), we shall give a justification of G.N. Hounsfield’s idea applying a stochastic approach.

In the system Equation (17), we divide the right hand side into there components as Equation (25);

$$\mathit{s}={\mathit{s}}^0+{\mathit{s}}_{\epsilon }^{app}+{\mathit{s}}_{\epsilon }^{obs},$$

(25)

where ${\mathit{s}}^0:=A{\mathit{x}}^0$ with Equation (26)

$${\mathit{x}}^0=(x_1^0,\cdots ,x_n^0),x_j^0:=\frac{1}{|c_j|}{\int }_{c_j}f(x,y)dxdy,$$

(26)

${\mathit{s}}_{\epsilon }^{app}:=({\left(s_{\epsilon }^{app}\right)}_1,\cdots ,{\left(s_{\epsilon }^{app}\right)}_m)$ with ${\left(s_{\epsilon }^{app}\right)}_j={\int }_{l_j}f(x,y)dl-s_j^0$ and ${\mathit{s}}_{\epsilon }^{obs}={\mathit{s}}_X:=(s_1^X,\cdots ,s_m^X)$ with $s_j^X,j=1,\cdots ,m,$ being random variables whose probability density functions are supported by $[-{\epsilon }_j,{\epsilon }_j]$ for small ${\epsilon }_j>0$ and whose expectation value $E(s_j^X)$ is 0. We have assumed that the $s_j^X$’s are the white noises, which is a very natural assumption. In this setting, ${\mathit{x}}^0$ is a theoretical solution to the system $A{\mathit{x}}^0={\mathit{s}}^0$ assigning the integral mean value of the function $x_j$ in each pixel $c_j$, that is, we take the pixel function as Equation (27)

$$g_n^0(x,y):=\sum _{j=1}^n{x}_j^0{\chi }_{c_j}(x,y).$$

(27)

The least square solution ${\mathit{x}}^{app}$ of the system $A\mathit{x}={\mathit{s}}_{\epsilon }^{app}$ gives the effect of the approximation of the function f by the pixel function Equation (27) to the least square solution of the system $A\mathit{x}=\mathit{s}$. The least square solution ${\mathit{x}}^{obs}$ of the system $A\mathit{x}={\mathit{s}}_{\epsilon }^{obs}$ give the effect of the observation errors, numerical rounding errors by computers and so on. We note that this assumption is natural since we can a priori estimate the observation errors, the size of the support ${\epsilon }_j$ of $s_j^X$ each for each j, and so on. In the same way as Equation (27), we define the pixel function determined by ${\mathit{x}}^{app}={}^t(x_1^{app},\cdots ,x_n^{app})$ and by ${\mathit{x}}^{obs}={}^t(x_1^{obs},\cdots ,x_n^{obs})$, respectively as Equation (28)

$$\begin{array}{c}\hfill g_n^{app}(x,y):=\sum _{j=1}^n{x}_j^{app}{\chi }_{c_j}(x,y),\\ \hfill g_n^{obs}(x,y):=\sum _{j=1}^n{x}_j^{obs}{\chi }_{c_j}(x,y).\end{array}$$

(28)

Let us fix the number n of the pixels to be large. By virtue of Theorem 2, the least square solution ${\mathit{x}}^{obs}$ of the system $A\mathit{x}={\mathit{s}}_{\epsilon }^{obs}$ is obtained by solving the system Equation (24) having the unique solution, which yields that each component of Equation (29)

$${\mathit{x}}^{obs}={\left({}^{t}AA\right)}^{-1}{}^{t}A\mathit{s}$$

(29)

can be regarded as the arithmetical mean of m radon variables whose expectations are 0. Therefore, the (strong) law of large numbers gives that Equation (30)

$$\underset{m\to \infty }{lim}{g}_n^{obs}(x,y)=0,a.s.,$$

(30)

which proves the following.

Theorem 2.

We know by experience that even for a fixed number n, the larger the number m is, the better resolution we can obtain. The Equation (30) proves that our knowledge by experience is stochastic theoretically right.

This theorem stochastically proves that the idea by G.N Hounsfield is theoretically right, in the sense that his solution approximates the unique least square solution to Equation (24) if the number n of the pixels are sufficiently large, we have enough number $m>n$ of the X-ray data and the errors are sufficiently small, which can be naturally assumed in usual practical applications.

3.4. A Proposal to Develop USCT

In this subsection, based on Property 1 and Hounsfield’s idea, we propose an idea for the development of USCT. The basis of our idea is that if the ultrasound rectilinearly propagates then Problems 1 and 2 are equivalent.

Theorem 3.

Assume the same assumption as Problem 1. If we assume that the ultrasound rectilinearly propagates then the fastest route between the two boundary points $\alpha ,\beta \in \partial \mathsf{\Omega }$ becomes the segment $l(\alpha ,\beta )$ as Equation (31)

$$\left\{\right(\alpha ,\beta )=\{t\alpha +(1-t)\beta |0<t<1\}$$

(31)

between α and β. In this case, the integral Equation (9)> the line integral along the segment between the two boundary points $\alpha ,\beta \in \partial \mathsf{\Omega }$. Therefore, if the ultrasound rectilinearly propagates then Problems 1 and 2 are equivalent by identifying $f\left(x\right)=1/C\left(x\right)$ for $x\in {\mathbb{R}}^2$ or ${\mathbb{R}}^3$.

By this theorem, if we assume the ultrasound rectilinearly propagates, then applying the same algorithm as X-ray CT, then we obtain some approximation for the section of the reciprocal velocity $1/C\left(x\right)$ distribution.

In this subsection, we shall propose more precise three dimensional reiterating reconstruction algorithm for $1/C\left(x\right)$ in USCT, where we shall apply Theorem 3 as the first step.

Method 1.

For USCT, we propose the following reiterating procedure.

(i)

Divide the object into boxels.

(ii)

Assuming that the ultrasound rectilinearly propagates and apply G.N. Hounsfield’s idea to obtain the boxel function $g_1$ as the least square solution to Equation (17), which approximates $1/C\left(x\right)$, the reciprocal of the sonic velocity $C\left(x\right)$. In this procedure, the elements $a_{ij}$ of the matrix in Equation (17) is given by the length of intersection between the linear (assumed) orbit of the i-th ultrasound and j-th boxel and $s_i$ is the travel time of the i-th ultrasound.

(iii)

Taking the boxel function $g_1$ obtained in the procedure (ii) as the first image. By the first image, we can calculate the new orbit of the ultrasound. Then we change the elements $a_{ij}$ of the matrix by the length of intersection between the new orbit of the i-th ultrasound and j-th boxel. Giving the least square solution $g_2$ to Equation (17) with new matrix elements gives better approximation.

(iv)

Taking the boxel function $g_2$ as the second image and repeat the same procedure as procedure (iii) and obtain a better approximation boxel function $g_2$.

(v)

Repeat the procedure (iv) until $sup|g_{k+1}-g_k|$ is sufficiently small.

Let us explain this method. In the first step (Method 1 (i), (ii)), we assume that the ultrasound rectilinearly propagates and apply the same algorithm as G.N. Hounsfield to obtain the boxel function $g_1$, which give the first image for our reiterating algorithm. In this step, we assume that the ultrasound rectilinearly propagates in the whole structure as well as in each boxel. Confer Figure 3 for this image where we give a two-dimensional picture for simplicity.

In the second step, having known the distribution function $g_1$ of the velocity, we can calculate the fastest route between the two boundary points. In this case, the route where the ultrasound propergates is not a line, in general, but the ultrasound rectilinearly propagates in each boxel (confer Figure 4). Therefore, the fastest route becomes a polygonal line and each component $a_{ij}$ in the overdetermined system Equation (17) (or Equation (18)) is the intersection of the fastest route of the i-th ultrasound and the j-th boxel $c_j$.

By changing the coefficients $a_{ij}$’s in the overdetermined system Equation (17) (or Equation (18)), we obtain the better approximation function $g_2$ of $1/C\left(x\right)$ than $g_1$, as the least square solution to the system Equation (17) (or Equation (18)).

In the third step, we calculate the new fastest routes based on $g_2$ and repeat the same procedure as the second step.

Reiterating the same procedure until $sup|g_{k+1}-g_k|$ is sufficiently small then $g_{k+1}$ is a sufficiently approximate boxel function of $1/C\left(x\right)$.

Therefore, in order to practicalize Method 1, it is very important to solve the following problem.

Problem 5.

Assume that the distribution $C\left(x\right)$ of the velocity, $x\in \mathsf{\Omega }\subset {\mathbb{R}}^3$ is known. For any $\alpha ,\beta \in \partial \mathsf{\Omega }$, determine the fastest route γ from α to β though Ω, that is, as Equation (32)

$$\underset{{\gamma }_{\alpha ,\beta }}{min}{\int }_{{\gamma }_{\alpha ,\beta }}1/C\left(x\right)d\gamma ={\int }_{\gamma }1/C\left(x\right)d\gamma .$$

(32)

This problem is challenging and interesting. The complete solution to Problem 5 seems to be difficult, although it is sufficient to solve this problem by assuming that the distribution $C\left(x\right)$ of the velocity is a boxel function. Furthermore, giving an algorithm to find the fastest route for a boxel function $C\left(x\right)$ is also practically sufficient.

The following algorithm is also sufficient as an alternative of the solution to Problem for boxel functions.

Method 2.

As for an alternative for the solution Problem for boxel functions, we propose the following procedure. Assume that we have a boxel function as Equation (33)

$$g\left(x\right)=\sum _{j=1}^n{x}_j{\chi }_{c_j}\left(x\right),$$

(33)

is obtained for $x\in {\mathbb{R}}^3$ (or $x\in {\mathbb{R}}^2$).

(i)

Find the boxes $c_j$’s where $x_j$’s are very large, which implies that the sonic velocity $1/x_j$ in the boxel $c_j$ is so small that the ultrasound would not go through it in practice.

(ii)

Find a box $c_k$ near $c_j$ whose value $x_k$ is larger than $x_j$.

(iii)

Change the route ${\gamma }_i$ of the i-th ultrasound by the new one ${\gamma }_i$’ where the ultrasound goes through $c_k$’s instead of $c_j$’s. We also change $a_{ij}$ to $a_{ij}$’ in accordance with the change of the routes.

(iv)

If ${\sum }_i{a}_{ij}{c}_j<{\sum }_i{a}_{ij}^{\prime }{c}_j$ then we replace the i-th route ${\gamma }_i$ by the new one ${\gamma }_i$’.

By application of this method, the number of the reiterating times may be larger, but it converges to obtain a suitable approximate boxel solution $g\left(x\right)={\sum }_{j=1}^n{x}_j{\chi }_{c_j}\left(x\right)$ for $1/C\left(x\right)$.

4. Conclusions

In order to conclude this paper, let us summarize.

• We gave a short survey on new tomographic techniques without X-ray or MRI.
• We gave a theoretical justification of the idea by G.N. Hounsfield to practicalize the X-ray CT.
• Applying Property 1 and the idea by G.N. Hounsfield, we have proposed a new algorithm for USCT.

It is our main purpose to propose a new methodology for USCT (Method 1), where we have created a new research task to solve Problem 5. It is interesting to study Problem 5 itself as well as its derived problems mentioned at the end of the previous section.

It is also interesting to compare the new tomographic techniques introduced in the second section by some numerical simulation, however, for the time being, for each new tomographic technique, there is no consolidated canonical methodology yet developed. Therefore, this comparison can be our future task. It will be interesting research, even under the current situation, to try to compare the new tomographic techniques using numerical simulation.