Using the Sonification for Hardly Detectable Details in Medical Images.

The inverse sonification problem is investigated in this article in order to detect hardly capturing details in a medical image. The direct problem consists in converting the image data into sound signals by a transformation which involves three steps - data, acoustics parameters and sound representations. The inverse problem is reversing back the sound signals into image data. By using the known sonification operator, the inverse approach does not bring any gain in the sonified medical imaging. The replication of the image already known does not help the diagnosis and surgical operation. In order to bring gains in the medical imaging, a new sonification operator is advanced in this paper, by using the Burgers equation of sound propagation. The sonified medical imaging is useful in interpreting the medical imaging that, however powerful they may be, are never good enough to aid tumour surgery. The inverse approach is exercised on several medical images used to surgical operations.

Introduction

A significant effort has been devoted in recent years to improve the quality of medical images used to surgery[1]. The Roentgen’s discovery of X-rays (1895) was followed by the computed tomography, magnetic resonance imaging, nuclear imaging, and ultrasound-positioned medical imaging used for diagnoses and surgery.

To our knowledge, we are the first to apply the sonification theory to uncover hidden details in medical images, such as vessels, organs or tumors that cannot be directly seen with the eye. The beginnings of the sonification theory are dated in 1952 when Pollack evaluates auditory displays as a visualization tool by using the information theory[2,3]. The International Community for Auditory Display Conference organized by Kramer in 1992 has generated great interest for this multi-disciplinary theory, from the science and technology to the arts[4,5]. Licht pursues the history of the sound art in 2007 by highlighting ancient art such as Sonic Youth and contemporary art that led to challenging applications, including the works of Christian Marclay, LaMonte Young, Janet Cardiff, Rodney Graham and Laurie Anderson among many others[6].

The nano-guitar built at Cornell University from crystalline silicon no larger than a single human blood cell, invites bacteria inside a pacient to sing so that it can be easily detectable by a stethoscope[7]. The quantum whistle produces a nano-scale sound capable of detecting oscillations in the superfluid gases predicted by quantum theory[8].

The sonification theory allows new perspectives in the diagnosis of diseases, such as the Alzheimers’s dementia[9] and therapies in body movements, such as walking, twisting, lifting the arms or moving the legs[10]. The inverse problem of sonification, that is, the reversal of sound samples in new images is less studied so far, to the best of our knowledge. This is due to the fact that the known sonification operator does not bring any improvement of the medical image, because the theory behind it is the linear theory of sound motion. This article introduces a new sonification operator based on the nonlinear Burgers theory of the sound motion. The new operator has proved its ability to solve the inverse problem of sonification and to obtain essential gains in improving the medical image.

The paper is organized as follows: Section 2 is devoted to description of the direct problem of sonification. A description of the new sonification operator based on the Burgers equation of sound propagation is presented in Section 3. The methodology is described in Section 4, the applications in Section 5, while Section 6 contains the conclusions.

Direct Problem of Sonification

The direct problem of sonification, as known in the literature[11-14], is based on a sonification operator S0 to transform the image point data D into sound signals Y0 as S0: D → Y0, S0: x(t) → y0(t0, x(t), p0), where x(t) is a 1D string of point data, t is the data time, t0 is the sonification time, and p0 ⊆ P0 is the set of sonification parameters.

The parameters P0 = {k0, Δ0, , α0, β0, ϕ0, ε0, g0, γ0, H0} include k0, the factor of time compressor on the interval T0 = T/k0, Δ0 ≥ 0, the factor of dilation, , the reference frequency, α0, β0 ≥ 0, the pitch scaling factors, ϕ0 ≥ 1, the power distortion factor, ε0 ≥ 0, the amplitude threshold, g0, the gain function, γ0, the decay parameter and H0, the timbral control function.

The variables of the data domain are t, t, T. The signal x(t) can be divided into segments of different length that do not overlap, being expressed as a sequence x(n) of N = T × f samples at the rate f of T seconds duration. The time points t split in time the segments x(t). A possible division in M segments of x(t) isfor t0 = 0 and t = T. The duration of each segment is T = t − t.

Each segment x(t) is sonified as a single event longer or shorter than T

The general form for the sonified signal y0(t0) iswhere x(Δ0t0′) is the mean free segment, and x(t) is the trend signal at the starting point for pitch modulation. Parameter Δ0 gives the length of the event . If Δ0 = k0 the adjacent events do not overlap but they can overlap for Δ0 ≤ k0.

To introduce control of timbre, the operator H0 acts as the sine function, sowhere a(t0) is the modulation amplitude, f is the base frequency for the pitch range of sonification, and b(t0) is a pitch modulator. The amplitude modulator is defined aswhere φ0 is the amplitude modulator. A half-wave rectification is included for exceeding a threshold ε0 around the mean of the amplitude, a half-wave rectification is included

New Sonification Operator

We are looking for a new operator to replace (3) based on the Burgers equation of sound propagation. A digital image B is seen as a collection of N pixels. We suppose that B is subjected to external vibration force f(t) expressed as a sum of harmonic force F(t), and the generation sound force F(t). The force F(t) has the role to build the sonification operator. The response of B to f(t) is a new configuration b of all points P ∈ B at the time t. The vibration of B is described by Burgers equation[15]where x = (x1, x2, x3) is the position vector, v = (v1, v2, v3) is the acoustic velocity vector, τ = t − x/c0 is the retarded time, t is time, c0 is the velocity of sound motion in the linear approximation, b = (b1, b2, b3) are the dissipation coefficients, ρ0 is density of medium, β = (β1, β2, β3) is nonlinearity coefficients. Details on the pulse propagation in nonlinear 1D media can be found in[16-19].

Equation (7) admits the cnoidal solutions[20]. These solutions are localized waves that conserve their properties even after interaction among them, and then act somewhat like particles. This equation and other equations of the same kind (Schrödinger, Korteweg–de Vries equations etc.) have an infinite number of local conserved quantities, an infinite number of exact solutions expressed in terms of the Jacobi elliptic functions (cnoidal solutions) or the hyperbolic functions (solitons), and the simple formulae for nonlinear superposition of explicit solutions.

Given a known force F(t), we determine F(t) such that the acoustic power W radiated from B to be minimum. The W is written aswhere v is the velocity verifying (7) and p the acoustic pressure vector, A is the area of the image, and the subscript T represents the Hermitian transpose[21].

The solutions v, i = 1, 2, 3 of (7) are expressed aswhere , l is a finite number of degree of freedom of the cnoidal functions, 0 ≤ m ≤ 1 is the moduli of the Jacobean elliptic function, ω is frequency and the phase, k1, k2, k3 are components of the wave vector[20]. In the following, we stop to l = 2, and we will see that there are no sensible improvements in solutions for l > 2 The function F(t) is determined from

The unknown parameters , j = 1, 2, 3, are find by a genetic algorithm which minimizes the objective function ϒ(P) written with respect to residuals of (7) and (10)with

The genetic algorithm is running until it is reached a non-trivial minimizer, which will be a point at which (11) admits a global minimum.

The quality of results depends on the values of ϒ. The required precision is taken to be six places after the decimal point. The genetic parameters are assumed to be as follow: number of populations 200, ratio of reproduction 1.0, number of multi-point crossover 1, probability of mutation 0.5, and maximum number of generations 500.

Once determined the function F(t), the sonification operator S(D, t) is written aswhere D = {d1, d2, …, d}, d ∈ R is the point data domain of the original image. The data matrix D is obtained from B by applying the converter C1 as C1(B) = D. The converter C1 is defined in the next Section. Data D is arranged as a matrix with arbitrarily number of elements.

In (13), is the point data domain of the sonified image depending on time. If duration of the sonification procedure is T seconds, then the final sonified image is

This matrix D is shown in Fig. 1b. Each element of the matrix may contain color or nuances, interfaces or borders separating the colors and nuances, lines, curves and other objects (Fig. 1a).

Figure 1

(a) A fictitious image with elements which may contain borders separating the colors and nuances, line and curves lines. (b) The matrix D.

Equation (7) is solved for sharp interface continuity conditions for both displacements and stresses. The reflections by the edges of the grids are removed by the Dirichlet and Neumann boundary conditions and a strongly attenuative buffer. The reflection coefficient is[22]where j is the degree of approximation, and θ is the incidence angle. For more than one component of displacement, the Dirichlet and Neumann conditions alternate components at the boundaries. When more than one boundary is nonreflecting, more solutions are added to eliminate multiple reflections.

After sonification, the point data domain of the sonified image may contain small blurred areas with cavities and white dots, due to the inaccuracies of the original images. The convertor C2 is applied to in order to fill these bad zones with colors and objects by prolonging through continuity of the solutions (9) in the adjacent areas and points. The converter C2 is defined in the next Section.

Methodology

The problem to be investigated and solved in this paper can be formulated in three steps as:

Given a digital image B, a converter C1 is applied to B to construct the data matrix D as C1(B) = D;

The sonification operator S is apply to D to obtain the sonified image as . If duration of the sonification procedure is T seconds, then the sonified image is .

A converter C2 is applied to to obtain the final sonified image B as .

Explanation of each path follows:

Path 1. The picture B is uniform sampled in a grid with equal size boxes B = {b1, b2, …, b}, b ∈ R. The sampling must be sufficiently fine with special care devoted to eliminate (or at least to keep under control) the sources of numerical errors.

The converter C1 is applied to B to obtain the point data D = {d1, d2, …, d}, d ∈ R (Fig. 1b) aswhere M(F) is an alphabet map containing four filters F, i = 1, 2, 3, 4, i.e. F1 the color and nuances filter, F2 the interfaces and borders separating the colors and nuances filter, F3 the line and curved lines filter, and F4 other objects that appears in the image filter (Fig. 1a). The filters are controlled by a code J of minimizing the resolution loss and improving of the noise performance.

The filters have ability to notice fine image details be it color or line, with no connection to pixel count or pixel density. Scheme of the alphabet M(F) is presented in Fig. 2.

Figure 2

Scheme of the alphabet M(F).

Each filter contains a number of sub-filters sF, i = 1, 2, 3, 4, k = 1, 2, 3 …, K.

Filter F1 contains 39 sub-filters sF1, k = 1, 2, 3 …, 39, for 39 color nuances. Each sub-filter captures a single color. Each color has a code α, k = 1, 2, 3 …, 39.

Filter F2 contains sub-filters sF2, j = 1, 2, …, . Each interface has a code β, k = 1, 2, …, , obtained by interpolating of two codes α, k = 1, 2, 3 …, 39.

The filter F3 contains a routine for linear, polynomial or nonlinear curve fitting along with validation of fit tests. The line or curve is introduced by points and the routine determines the most suitable equation for it.

The filter F4 identifies other objects from B. The shape Γ of the object is defined as the image of the unit n-sphere S of equationthrough the affine transformationwhere r = r(ξ, ψ, ζ) are the components of rotation which transforms the coordinate axes into the principal axes of the sphere.

An inverse problem is applied to find the set of parameters (shape parameters) that define Γ, i.e. arbitrary center coordinates x, y, z, principal axes a, b, c, the principal directions defined by Euler angles ξ, ψ, ζ and the exponent n. The advantage of this model is the small number of parameters needed to represent a shape.

Each filter works across a different algorithm.

2. The sonification operator S(D, t) defined by (13) is applied to D to obtain the sonified image aswhere F(t) is the generation sound force determined from the condition of minimum acoustic power W. By setting , the function F(t) is determined aswith

In (16) is the point data domain of the sonified image, and t is the sonification time.

3. After sonification, may contain small blurred areas with cavities and white dots due to the inaccuracies of the original images. We term these areas as damaged zones.

The convertor C2 has the role to fill the damaged zones with color and geometric lines, through continuity of solutions (9) of adjacent areas and neighboring points.

The converter C2 is applied to to obtain the final sonified image B = {b1, b2, …, b}, b ∈ R, as wheremeasures the distance between the solution in a damaged point and the solutions in its neighboring points. In (23), δX is the solution in the damaged point, and δP is the prolonged solution through continuity of the solutions in neighboring points. determines the best solution in the damaged areas or points. Less spatial and color artifacts and better noise performance is assured by C2 compared with the existing schemes.

The filters are controlled by a code J of minimizing the resolution loss and improving of the noise performancewhere ε, ε, ε measure the relative errors for identification of colors, the lines and curves, and the moment of inertia of objects, defined as

In (24), α, α, α are the associated weights.

Applications

A sample of a fictive rat liver which exhibit changes in profile by severe loss of architecture and disturbances zones (between 10 and 50 μm) at the microscopic scale[23], is shown in Fig. 3a. The size of constituents is displayed in Fig. 3b.

Figure 3

(a) A fictive rat liver sample. (b) The constituent’s size.

The sonification operator (13) is exercised on fictive images of fibrotic rat liver samples inspired from an investigation of the effects of an extract of ginkgo biloba leaf against hepatic toxicity induced by methotrexate in rat[23,24]. The cross-sections of the rat liver are shown in Fig. 4.

Figure 4

Nine sectional slices of a rat liver sample. (b) Digital images of cross-sectional slices.

For the first two applications, the simulations were performed for digital images of 1.600 × 1.200 pixels (length and width). Sonification time is 6 sec., n = 3 for the unit n-sphere S, and the weights in the code J (α, α, α) = (0.4, 0.4, 0.2).

The properties for the rat liver are: density ρ0 = 1.05 g/cm3, dissipation coefficients b = (b1, b2, b3) = (0.2, 0.2, 0.2) kg sec/m2, nonlinear coefficients β = (β1, β2, β3) = (0.3, 0.3, 0.3) sec.

Figure 5 visualizes the new images obtained after sonification. By comparing these images to the original ones, some differences are highlighted in yellow in the last six images. Although there appears a tendency for the replication of images already known, our results show relevant details absent in the original images.

Figure 5

Cross-sectional slices of the sample after sonification. Yellow circles indicates the differences with the original images.

For the next exercise, we consider the work of Salameh[25] which studies the detection of nonalcoholic steatohepatitis in the fatty rat livers by magnetic resonance (MR). This study is useful in the early detection of fibrosis in the at livers[26-28]. Figure 6a shows the MR image of a liver rat with strong hepatocellular damages. Some details are purposefully hidden (red circles in Fig. 6b). We see that the inverse sonification operator recovered all initially hidden details (Fig. 6c).

Figure 6

(a) The MR image of a liver rat. (b) Initially hidden details - red circles. (c) The hidden details were recovered by the sonification technique.

Another application is related to the hepatic arterial chemotherapy[29-32]. The knowledge of hepatic and biliary vascular maps is absolutely necessary for planning the surgical operation. A catheter must be inserted inside the gastroduodenal artery (GDA) to distribute the chemotherapy. A possible location of the hepatic arterial infusion catheter was discussed in[30] and shown in Fig. 7.

Figure 7

A possible location of the hepatic arterial infusion catheter[30].

Figure 8a shows the CT image of the hepatic artery (CHA- common hepatic artery, LHA - left hepatic artery, RHA - right hepatic artery, SA - splenic artery, Seg IV HA - segment IV hepatic artery)[30]. Figure 8b shows the CT image of the left hepatic artery[30].

Figure 8

(a) CT image of the hepatic artery (CHA - common hepatic artery, LHA- left hepatic artery, RHA - right hepatic artery, SA - splenic artery, Seg IV HA - segment IV hepatic artery). (b) CT image of the left hepatic artery; (c) the image before sonification; (d) final sonified image.

Our intension is to obtain an image by sonification of the image shown in Fig. 8c.

For last applications, the simulations were performed for digital images of 1.600 × 1.200 pixels (length and width). Sonification time is 7 sec., n = 3 for the unit n-sphere S, and the weights in the code J (α, α, α) = (0.4, 0.4, 0.2).

The properties for the rat liver are: density ρ0 = 1.12 g/cm3, dissipation coefficients b = (b1, b2, b3) = (0.2, 0.2, 0.2) kg sec/m2, nonlinear coefficients β = (β1, β2, β3) = (0.3, 0.3, 0.3) sec.

Final sonified image is shown in Fig. 8d. We see that it is identically to the one shown in Fig. 7.

The last application considers the case of a tumor (pink color) located near the portal tree of the vascular territory (Fig. 9a)[33,34]. The vascular territory (1) and the vessel branches in the vicinity of tumor (2) are shown in Fig. 9b.

Figure 9

(a) The tumor location. (b) Vascular territory (1) and the vessel branches in the vicinity of the tumor (2).

After sonification three images are obtained in Fig. 10 for the frontal, caudal and cranial views. New details on the tumor and surrounding areas are obtained and aditionally, the shape and size of the tumor is better visualized.

Figure 10

Sonified images in the vicinity of the tumor after sonification.

Conclusion

The sonification operator proposed in this paper is converting the digital data field into sounds by using the nonlinear Burgers sound equation which is solved by cnoidal method. By inverting the sound into image, the result highlights hidden details in the image seen by the sound and not seen by the eyes. To show the efficiency of the sonification procedure and to verify the correctness of the results, we intentionally hide some details in the images before the sonification. The sonification operator has a positive effect in improving the accuracy in the medical imaging by capturing and detecting hard-to-see details in original images.

The advantages of the sonification refer to: (1) there is no distinction between the data time and the sonification time; (2) reflects indirectly the relationship between the propagation of sound through tissues and the structure of these tissues; (3) discovers new hard-to-find details in original medical images; (4) demonstrates the usefulness of image-sound and sound-image conversions in exploiting the medical imaging in diagnosis and surgery.