VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ BRNO UNIVERSITY OF TECHNOLOGY FAKULTA ELEKTROTECHNIKY A KOMUNIKAČNÍCH TECHNOLOGIÍ ÚSTAV BIOMEDICÍNSKÉHO INŽENÝRSTVÍ FACULTY OF ELECTRICAL ENGINEERING AND COMMUNICATION DEPARTMENT OF BIOMEDICAL ENGINEERING CALIBRATION OF AN ULTRASONIC TRANSMISSIVE COMPUTED TOMOGRAPHY SYSTEM KALIBRACE ULTRAZVUKOVÉHO PRŮZVUČNÉHO SYSTÉMU VÝPOČETNÍ TOMOGRAFIE DOKTORSKÁ PRÁCE DOCTORAL THESIS AUTOR PRÁCE Ing. ADAM FILIPÍK AUTHOR VEDOUCÍ PRÁCE prof. Ing. JIŘÍ JAN, CSc. SUPERVISOR BRNO 2009 Keywords Transmissive ultrasonic computed tomography, ultrasound attenuation imaging, sensi- tivity calibration, geometrical calibration Klíčová slova Průzvučná ultrazvuková počítačová tomografie, zobrazování ultrazvukového útlumu, kalibrace citlivostí, geometrická kalibrace ii Acknowledgements This dissertation thesis would have not been possible without the support of several people. First of all, I want to thank my academic advisor Prof. Jiří Jan for giving me direction throughout the whole doctoral study. His observations made me rethink and reformulate large parts of the dissertation work making it more clear and concise. A special thank goes to all members of the Karlsruhe USCT team, namely Dr. Rainer Stotzka, Dr. Nicole Ruiter, and Michael Zapf, who gave me excellent support during the months I spent in the Karlsruhe Forschungszentrum. Further I would like to thank all members of the DAR team, who formed a very good basis for exchanging ideas. Especially Radovan Jiřík PhD gave me a lot of advice in the initial stages of work and also formed the route I would take in this dissertation. Igor Peterlík PhD and Dušan Hemzal PhD I want to thank for helping me mathematically formulate some of the most important parts of the dissertation. Many thanks also go to my fellow doctoral students making the studies enjoyable and unforgettable. And last but most importantly I want to thank my wife and my parents for giving me support throughout my whole studies. Thank you. iii Frame of the dissertation The presented dissertation project has been running from September 2003 to August 2009 at the Department of Biomedical Engineering, Brno University of Technology, under the academic supervision of Prof. Jiří Jan. The work direction was determined by the long-term orientation of research in the department towards ultrasonic related problems and especially by the research cooperation with the Institute of Data Processing and Electronics, Forschungszen- trum Karlsruhe (German National Research Center). A substantial part of the experimental research was done in the Karlsruhe research fa- cility, where the author joined the development team several times for a total of 4 months. The team in Karlsruhe is developing a prototype of a novel 3D ultrasound tomography system aimed at breast cancer screening. The cooperation and the research visit were great sources of inspiration, exchange of ideas, and most importantly experimental data, on which newly de- veloped methods could be tested. One of the main problems in the Karlsruhe USCT prototype was the absence of any ca- libration. The thousands of transducers used in the system have deviations in sensitivity, direc- tivity, and frequency response. Moreover, these parameters change over time as the transduc- ers age. Finally (and most importantly), the mechanical positioning of the transducer elements is not precise. All these aspects greatly affect the overall quality of the reconstructed images. The problem of calibration of a USCT system was chosen as the main topic for this disserta- tion. The dissertation work and the stay in the Karlsruhe research facility was supported by:  National research center DAR (Data, Algorithms, and Decision making), sponsored by the Ministry of Education, Czech Republic, project no. 1M6798555601  Bilateral Czech-German agreement AVČR - DAAD (grants no. D-CZ 22/05-06 and D12-CZ9/07-08 and the grants of the German Academic Exchange Service )  Research frame of the FEEC BUT (Ministry of Education, Czech Republic), project no. MS 0021630513 (MS 1850023) iv Used Acronyms  DOP – dilution of precision (described in detail in chapter 2.4.1.4)  FZK – Forschungszentrum Karlsruhe  GPS – global positioning system  ITE – individual transducer element, an alternative approach to the novel geometrical cali- bration technique (described in chapter 6.1.1)  TAS – transducer array system, a module consisting of several ultrasound transducer ele- ments (chapter 2.2.2.2) or also an alternative approach to the novel geometrical calibration technique (described in chapter 6.1.2)  TOA – time-of-arrival  TOF – time-of-flight  USCT – ultrasound computer tomography  XCT – X-ray computer tomography 1 Contents ACKNOWLEDGEMENTS ........................................................................................................................ II FRAME OF THE DISSERTATION ............................................................................................................ III USED ACRONYMS ............................................................................................................................... IV 1. INTRODUCTION .......................................................................................................................... 5 1.1 COMPUTED TOMOGRAPHY SYSTEMS ...................................................................................................... 5 1.2 ULTRASOUND COMPUTED TOMOGRAPHY ............................................................................................... 6 1.2.1 Example of a USCT system .................................................................................................. 6 1.2.2 History of USCT ................................................................................................................... 7 1.3 IMAGE RECONSTRUCTION PRINCIPLES .................................................................................................... 8 1.3.1 Filtered backprojection ....................................................................................................... 8 1.3.2 Diffraction tomography .................................................................................................... 10 1.3.3 Reflectivity imaging........................................................................................................... 11 1.3.4 Algebraic Reconstruction Techniques ............................................................................... 12 1.4 THE NEED FOR CALIBRATION .............................................................................................................. 12 2 STATE OF THE ART IN USCT ....................................................................................................... 14 2.1 ULTRASOUND COMPUTED TOMOGRAPHY SYSTEMS ................................................................................. 14 2.1.1 High-Resolution Ultrasonic Transmission Tomography System ........................................ 14 2.1.2 The Ring Transducer System for Medical Ultrasound Research ........................................ 14 2.1.3 A diffraction tomography system ..................................................................................... 15 2.1.4 Computerized Ultrasound Risk Evaluation (CURE) program ............................................. 16 2.1.5 Karlsruhe USCT system ...................................................................................................... 17 2.2 KARLSRUHE USCT SYSTEM AND PREVIOUS RESULTS ................................................................................ 17 2.2.1 The Karlsruhe 2D USCT system.......................................................................................... 17 2.2.2 The Karlsruhe 3D USCT system.......................................................................................... 18 2.2.2.1 System architecture .................................................................................................................... 18 2.2.2.2 Transducer array systems - TAS .................................................................................................. 19 2.2.3 FZK Reconstruction method (the sum-and-delay algorithm) ............................................ 20 2.3 PUBLISHED SENSITIVITY CALIBRATION METHODS ..................................................................................... 22 2.4 KNOWN POSITION CALIBRATION METHODS ........................................................................................... 23 2.4.1 Global Positioning System (GPS) ....................................................................................... 23 2.4.1.1 Pseudorange equations............................................................................................................... 24 2.4.1.2 Linearization ................................................................................................................................ 24 2.4.1.3 Solving the system of equations ................................................................................................. 25 2.4.1.4 Dilution of Precision (DOP) .......................................................................................................... 27 2.4.1.5 The use of GPS principles in USCT ............................................................................................... 28 2.4.2 A published GPS modification for an ultrasonic systems .................................................. 28 2.4.2.1 The calibration method ............................................................................................................... 29 2.4.3 The RTS calibration approach ........................................................................................... 31 2.4.4 Multidimensional scaling .................................................................................................. 31 2.4.5 Localization methods in wireless networks ....................................................................... 33 3 AIMS OF THE DISSERTATION .................................................................................................... 34 4 A CONTRIBUTION TO ATTENUATION IMAGE RECONSTRUCTION .............................................. 35 4.1 ATTENUATION COEFFICIENT CALCULATION ............................................................................................ 35 4.2 THE STRAIGHT BEAM RECONSTRUCTION ............................................................................................... 36 4.3 THE REFLECTED BEAM RECONSTRUCTION .............................................................................................. 38 4.4 RECONSTRUCTION RESULTS ............................................................................................................... 39 5 NOVEL TRANSDUCER SENSITIVITY CALIBRATION METHOD ....................................................... 41 5.1 THE 2D SENSITIVITY CALIBRATION METHOD .......................................................................................... 41 2 5.1.1 The measured signal model ............................................................................................. 41 5.1.2 The equation system ........................................................................................................ 42 5.1.3 Stabilizing the system ....................................................................................................... 42 5.1.4 Sensitivity calibration results based on experimental data .............................................. 44 5.1.4.1 Verification via hydrophone measurement ................................................................................ 45 5.1.4.2 Comparison with wave-equation based simulation ................................................................... 46 5.2 THEORETICAL EXTENSION TO 3D ........................................................................................................ 47 6 NOVEL TRANSDUCER POSITION CALIBRATION METHOD .......................................................... 49 6.1 FORMULATION OF THE POSITION CALIBRATION METHOD ......................................................................... 49 6.1.1 The ITE (individual transducer element) approach........................................................... 50 6.1.2 The TAS (transducer array system) approach .................................................................. 53 6.2 ANCHORING .................................................................................................................................. 56 6.2.1 Position anchoring............................................................................................................ 57 6.2.1.1 Dilution of precision (DOP) ......................................................................................................... 57 6.2.1.2 Minimizing the DOP for a 2D USCT ............................................................................................. 58 6.2.1.3 Minimizing the DOP for a 3D USCT ............................................................................................. 59 6.2.2 Time-delay anchoring ....................................................................................................... 63 6.2.3 The rank deficiency test.................................................................................................... 64 6.3 NUMERICAL ANALYSIS AND TESTING ................................................................................................... 64 6.3.1 Convergence analysis ....................................................................................................... 64 6.3.2 Condition number analysis ............................................................................................... 66 6.3.3 Noise effect analysis ......................................................................................................... 67 6.4 ALTERNATIVE SOLUTIONS TO THE CALIBRATION PROBLEM ....................................................................... 69 6.4.1 The Levenberg-Marquardt algorithm .............................................................................. 69 6.4.2 The weighted least squares approach.............................................................................. 70 6.4.3 Solution based on a common time-delay parameter ....................................................... 70 6.4.4 Solution based on angle-dependent common time-delay ................................................ 71 6.4.5 Theoretical extension to a thread-phantom-based calibration ....................................... 71 6.4.5.1 2D thread reflections .................................................................................................................. 71 6.4.5.2 3D thread reflections .................................................................................................................. 72 6.4.5.3 Expressing the TOA equations .................................................................................................... 74 6.4.5.4 Phantom description .................................................................................................................. 75 6.5 POSITION CALIBRATION OF THE KARLSRUHE 3D USCT SYSTEM ................................................................ 76 6.5.1 Extracting the TOA from A-scans ..................................................................................... 76 6.5.1.1 Pulse detection ........................................................................................................................... 77 6.5.1.2 Pulse selection ............................................................................................................................ 79 6.5.2 Analysis of measurement errors ....................................................................................... 83 6.5.2.1 Speed of sound errors ................................................................................................................ 83 6.5.2.1.1 Thermometer errors .......................................................................................................... 84 6.5.2.1.2 Temperature variations ...................................................................................................... 84 6.5.2.1.3 Error of the speed of sound ............................................................................................... 85 6.5.2.2 Thermal expansion errors .......................................................................................................... 86 6.5.3 Calibration results based on experimental data .............................................................. 87 6.5.3.1 The ITE and the TAS approaches ................................................................................................ 87 6.5.3.2 Different TOA selection methods ............................................................................................... 88 6.5.3.3 Alternative calibration solutions ................................................................................................ 88 7 CONCLUSION ............................................................................................................................ 93 7.1 DISCUSSION OF THE ACHIEVED RESULTS ............................................................................................... 93 7.1.1 Attenuation image reconstruction ................................................................................... 93 7.1.2 Sensitivity calibration ....................................................................................................... 93 7.1.3 Geometrical calibration .................................................................................................... 94 7.2 SUGGESTIONS FOR FURTHER RESEARCH ............................................................................................... 95 7.2.1 Ethanol solution measurement medium .......................................................................... 95 7.2.2 3D extension to sensitivity calibration ............................................................................. 95 7.2.3 Thread phantom calibration ............................................................................................ 95 7.2.4 LSQR solution to the system of equations ........................................................................ 96 7.2.5 Different pulse shapes, Pulse trains, Chirps ..................................................................... 96 3 7.2.6 Single- and double-differences .......................................................................................... 96 8 REFERENCES ............................................................................................................................. 97 9 LIST OF AUTHOR’S PUBLICATIONS .......................................................................................... 101 4 5 1. Introduction This dissertation is centered on a medical imaging modality – the ultrasonic computed tomography (USCT) – and algorithms which improve the resulting image quality, namely the calibration of a USCT device. The USCT is in principle capable of producing quantitative 3D image volumes with high resolution and tissue contrast. Although the idea of USCT has been around for more than three decades, the lack of computing power did not allow the scientists and engineers to make full use of its potential. With computers that are more and more powerful, the idea of USCT came back to life and at the moment several teams around the world are experimenting on prototypes of USCT scanners. The USCT is primarily aimed at breast cancer diagnosis. Breast cancer is the most common cancer among women worldwide. In 2000, the last year for which global data exists, over one million cases have been diagnosed around the world and some 400.000 women have died, representing 1.6 per cent of all female deaths. Early detection of breast cancer is vital since early detection has repeatedly been shown to improve the chance of survival [74]. X-ray mammography, magnetic resonance, and the conventional ultrasound are estab- lished methods for breast cancer diagnosis. Currently the most common modality used for breast cancer diagnosis is X-ray, but is often supplemented with an ultrasound examination, which in many cases leads to additional information about e.g. cysts and fibroadenomas. Be- sides the fact that ultrasound is much less expensive than X-ray mammography, it can be ap- plied more frequently. There is still an ongoing debate about the limits on the energy which the human tissue should be exposed to by ultrasonic waves (especially the thermal damage to the nervous system is a potential risk [6]). Ultrasound waves are however not ionizing and can therefore be applied on a regular basis. The disadvantage of the conventional echo ultrasound methods is mainly poor resolu- tion and therefore the inability to reliably detect microcalcifications. Both the resolution and contrast depend on the distance of the ultrasonic probe from the region of interest inside the breast. As the medical doctor tries to get as close to the desired area, the breast gets deformed, and therefore exact measurements of the tissue structure, e.g. the tumor size are hardly possi- ble. The image contents and their quality are highly dependable on the expertise of the ex- aminer and are hardly reproducible. In the following subchapters of the introduction, the principles USCT will be ex- plained. First a brief overview of the kinds of existing computed tomography systems (used in medical imaging) is given in chapter 1.1. Chapter 1.2 gives explanations on the principles of the USCT and some ultrasound-specific aspects of this imaging modality. The general prin- ciples of image reconstruction methods are explained in chapter 1.3. Finally, chapter 1.4 short- ly describes why there is a need to calibrate USCT systems and what can be calibrated. 1.1 Computed tomography systems Tomographic systems, are systems providing images by sections. Computed tomogra- phy systems are those, which feed the collected data to a tomographic reconstruction software yielding the tomographic images after being processed by a computer. The modern systems usually reconstruct a whole 3D volume and then offer the examiners to view a 2D slice of the volume at any angle. The computed tomography systems currently used for medical imaging can be classi- fied into the following categories according to the imaging medium or the basic imaging prin- ciple [29]: 6  X-ray computed tomography (XCT)  Nuclear imaging (SPECT and PET)  Ultrasonic imaging (USG)  Magnetic resonance imaging (MRI)  Electrical impedance tomography (IT) Although the conventional B-mode ultrasonography is usually not considered a com- puted tomography system, because in principle it does not need a computer to yield sectional images, the generation of images in all modern systems is assisted by a computer. This is es- pecially true, when considering the new three-dimensional ultrasonic systems or freehand ul- trasound [7], [29]. The ultrasonic computed tomography (USCT) is a completely new approach combin- ing the conventional ultrasonic imaging with the principles of image reconstruction used in XCT. 1.2 Ultrasound computed tomography The ultrasonic computed tomography is an imaging modality which combines the phe- nomenon of ultrasound and some image reconstruction principles developed for other tomo- graphic systems. The basic objective of a USCT system is the same as for other tomographic equipment: to obtain scans of all possible directions around the object and then reconstruct the entire volume. In comparison with the conventional B-mode ultrasound, the USCT receiving transducer elements need not necessarily be in the same place as the emitting elements. There- fore these systems are sometimes referred to as transmissive ultrasound systems. 1.2.1 Example of a USCT system In practice, a USCT system can resemble the following: thousands of ultrasonic trans- ducers attached to the inner side of a circular frame encompass the imaged object. Since the ultrasonic waves are not transmissible in the air (at the frequencies used in diagnostic medi- cine), water has to be used as coupling medium to create an acoustic link. Figure 1-1: A schematic drawing of the Karlsruhe 3D USCT prototype. (Taken from [69]) Scanning is done in the following way: one emitter at a time broadcasts a broadband ultrasonic pulse wave and all receivers record the direct, reflected and scattered ultrasonic waves. This procedure is repeated for each emitter separately creating a fan projection each time (Figure 1-2). Several ultrasonic properties of the scanned object can be estimated from each of the acquired signal can. These are: attenuation, reflectivity and the speed of sound. Here we can find a difference from the X-ray tomography in which it is only possible to calculate the atten- 7 uation of the ray intensity as it passes through the human tissue. An XCT detector integrates the radiation over a small period of time resulting in one value of intensity. The USCT receiv- ers however, can record longer signals containing a series of immediate pressure values. Figure 1-2: Architecture of the 2D USCT system built in Karlsruhe. A ring of ultrasonic transducers surround the examined object. One transducer transmits an ultrasonic pulse while all other transducers receive simulta- neously. (Taken from [65]) Thus, USCT offers significantly more information about the scanned object. At the same time, this approach poses great problems as the volume of the recorded data dramatically increases and one has to deal with increased demands on storage space and processing speed. 1.2.2 History of USCT The idea of ultrasound tomography is not new. The earliest attempts were made in the late seventies and early eighties of the last century. One of the first works in this area [21] dealt with only reconstructing the absorption of the tissue, the same property (in terms of sig- nal processing) as the one reconstructed in the XCT systems. Soon, however, it was clear that the USCT systems have a much greater potential. The speed of the ultrasound waves is rela- tively small (compared to the speed of electromagnetic waves - the speed of light) and easily measured. Because it varies in different kinds of tissues, the local velocity is an additional pa- rameter which can be reconstructed [22]. Norton and Linzer [52] showed that it is possible to reconstruct the “reflectivity” (like the conventional B-mode systems). Although reflectivity is not a precisely defined property of the tissue and corresponds to a combination of scattering, refraction and reflection, the images provide information about the structure of the tissue and have a high diagnostic value. This work also developed the idea of using unfocused transducers, which emit into the whole scanned volume. In previous experiments only transducers radiating a narrow focused ultra- sonic beam were used. This change corresponds to the transition of XCT systems from the first to the third generation. Making one projection of the entire area at once, rather than se- quentially, shortens the minimum period of capturing the entire scene from several hours to several seconds. In spite of the limited technology available to these investigators, they showed promis- ing results. For example, Greenleaf et al. achieved a sensitivity of 100% for palpable lesions with USCT for a small sample population [23]. In a larger study (n=78), Schreiman et al. showed that a computer-aided diagnosis using USCT had a sensitivity of 82.5% for the diag- nosis of a malignancy [62]. One of the biggest problems encountered was the enormous amounts of data flows, which the technology was not able to cope with. Especially the researchers were not able to record enough data (at least not quickly enough) [34]. And so after the first years of significant 8 development, USCT gradually retreated from the practical implementation (the companies Philips and Siemens closed their programs of developing commercial systems), and further work continued primarily on the theoretical level [4]. In the last decade or two, however, the researchers have again slowly gained interest and the topic of USCT has been reopened. The most recent approaches are documented in chapters 2.1 and 2.2. 1.3 Image reconstruction principles The first approaches to reconstruct images of the USCT led naturally in the footsteps of the already developed methods used in XCT. These methods assume the radiation of signals along straight thin lines - rays. The phenomena such as reflection, refraction and scattering are not taken into account. In such cases the filtered backprojection can be used as an effective method to reconstruct attenuation and velocity images (chapter 1.3.1). A different approach the diffraction tomography allows the “rays” to be slightly bent (1.3.2). Finally, reflectivity imaging takes advantage of the whole length of the recorded ultrasonic signals (chapter 1.3.3). 1.3.1 Filtered backprojection The simplest (and a very fast) method of reconstruction is called the backprojection method. The method assumes that the imaging system collects data organized in the so-called projections 𝑝𝜃0 𝜏 = 𝑓(𝑥,𝑦)𝑑𝑟 𝑅𝜏 . (1.1) Any projection 𝑝𝜃0 (𝜏) at an angle 𝜃0 consists of a set of ray-integrals (Figure 1-3). Each ray-integral is an individual measurement of a two-dimensional function 𝑓(𝑥,𝑦) (representing the distribution of a certain parameter to be imaged) integrated along a straight line. Figure 1-3: Ray integrals and a projection representation of an image. (Taken from [29]) If we have a continuous two-dimensional function (a so called sinogram) 9 𝑝 𝜏,𝜃 = 𝑓(𝑥,𝑦)𝑑𝑟 𝑅𝜏 ,𝜃 , 𝜏 ∈ −∞,∞ , 𝜃 ∈ (0,2𝜋) (1.2) both in 𝜃 and 𝜏, it is possible to invert it using the inverse Radon transform, and obtain the original function 𝑓(𝑥,𝑦) - the reconstructed image. Although the inverse Radon transform shows that, the original image can be reconstructed from its projections, it is very unstable in the presence of measurement noise. The backprojection algorithm is a practically usable and a stable alternative to the in- verse Radon transform. It follows a simple logic that each projection can be “smeared” back into the image in the direction originally used to for the projection. The resulting image is then the composition of all the smeared projections, where every point (x, y) is contributed to by a single value (the ray-integral) from each projection. 𝑏 𝑥, 𝑦 = 𝑝𝜃 𝜏 𝑑𝜃 = 𝜋 0 𝑝𝜃 𝑥 cos 𝜃 + 𝑦 sin𝜃 𝑑𝜃 𝜋 0 (1.3) The backprojection is not equivalent to the inverse Radon transform. Because the pro- jections are smeared over the whole reconstruction plain (𝑥,𝑦), an original image consisting of a single point will be reconstructed with non-zero values also outside of that point. A mod- ification of this algorithm was developed to account for this difference and is known as the filtered backprojection 𝑓 𝑥,𝑦 = 𝑞𝜃 𝜏 𝑑𝜃 = 𝜋 0 𝑞𝜃 𝑥 cos 𝜃 + 𝑦 sin𝜃 𝑑𝜃 𝜋 0 (1.4) where the modified projection 𝑞𝜃 𝜏 = 𝑝𝜃 ∗ 𝑕|(𝜏) is the convolution of the original projection with a ramp filter 𝑕 having a ramp-like frequency response |𝑤| (linearly growing with the fre- quency). The derivation the filtered backprojection is closely related to the Fourier slice theo- rem and is documented in numerous publications [7],[10],[12],[27],[28],[29],[35]. The backprojection algorithm can also be implemented in the frequency domain using the so called slice theorem. The slice theorem relates the spectrum of a two-dimensional func- tion to the spectrum of its parallel projection: ℑ 𝑝𝜃 𝜏 = ℑ2𝐷 𝑓 𝑥, 𝑦 𝜔 cos 𝜃 ,𝜔 sin𝜃 (1.5) That is, the spectrum of the projection 𝑝𝜃(𝜏) obtained from the image 𝑓 𝑥.𝑦 is equal to the central slice, at the angle 𝜃, of the two-dimensional image spectrum [29]. This theorem can be used to reconstruct an image from its projections in the following way. Each parallel projection is first transformed (by DFT) into its one-dimensional spectrum. The spectra are then filled according to the theorem (as central slices) into a two-dimensional function which becomes the two-dimensional spectrum of the image. Interpolation is needed to resample the 2D spectrum from polar to rectangular grid. Finally a 2D inverse DFT is ap- plied yielding the reconstructed image. The USCT devices usually do not provide data in the form of projections with parallel ray paths. In modern USCT systems, the transducers are usually placed on a circular frame. When the emitting element sends out a pulse wave it is recorded at once by all the receivers as fan-shaped equiangular projection. In order to be able to use the backprojection algorithm (which is only suited for the pa- rallel projection geometry) a procedure called rebinning is necessary to apply to the collected data. The main idea of rebinning is that each ray-integral from a fan projection 𝑟𝛽(𝛼) can also 10 be a part of a parallel projection 𝑝𝜃(𝜏).Thus the equiangular fan-projections can be rearranged into parallel projections and subsequently used in the backprojection algorithm to reconstruct the image. Attenuation imaging The filtered backprojection algorithm is for instance suitable to reconstruct images of local attenuation of the scanned tissue. Ultrasonic attenuation is derived from the recorded signals. There are three basic ways of determining the attenuation: the energy ratio method, the frequency shift method, and the method of log-spectra differences [14]. All three methods require a reference measurement, which can be done by making a scan with only water inside the tank. Then the imaged object is inserted into the system and the scan is repeated. The attenuation values are then computed using both scans. Each calcu- lated attenuation value represents the complete attenuation along the line from the emitter to the receiver - the ray-integral – and can be directly used as the input of the filtered backprojec- tion algorithm, yielding an image representing the distribution of local attenuation in the scanned object. Velocity imaging The distribution of local ultrasonic velocity can be reconstructed in a similar manner as the attenuation. The velocity can be simply calculated by taking the distance between the emit- ter and receiver and dividing it by the measured time-of-flight of the ultrasonic pulse. Al- though this method is the most straightforward, it is not very easy to determine the exact mo- ment of the arrival of the pulse. Rather than measuring the absolute velocities, one can obtain relative changes by making two scans (an empty one and one with the scanned object inside) [25]. No matter which method is used, the calculated speed of sound values (average veloci- ties along the ultrasonic “rays”) can be again used to reconstruct an image of the local veloci- ties inside the object using the filtered backprojection algorithm. 1.3.2 Diffraction tomography The above-mentioned reconstruction method assumes the radiation of the ultrasonic signals along straight thin rays. This assumption is not always valid, especially in the cases where the size of displayed objects is approximately same as the size of the ultrasound pulse wavelength. In such an environment, diffraction becomes a significant factor in the propaga- tion of the waves. The wave equations are then the more suitable means to describe the propa- gation of waves [35]. The solution of the wave equations is usually only possible under certain simplifica- tions. The two most cited are the Born and Rytov approximations (both simplify the equations by limiting the amount of diffraction to some small amount - about 10%) [35]. An important step in the reconstruction is the use of the Fourier diffraction theorem, which plays the same role in diffractive tomography as the slice theorem does in the conventional tomography. The theorem (1.6) says that the Fourier spectrum of a parallel projection (after illumi- nation of the object by a plane wave) is equal to the values of the two-dimensional spectrum of the image along a half-circle curve: 11 ℑ 𝑝𝜃 𝜏 = ℑ2𝐷 𝑓 𝑥,𝑦 𝜔 cos 𝜃 − 𝐾0 2 − 𝜔2 − 𝐾0 sin𝜃 ,𝜔 sin 𝜃 − 𝐾0 2 − 𝜔2 − 𝐾0 cos𝜃 (1.6) where 𝐾0 = 2𝜋 𝜆 , and 𝜆 is the wave length. The radius of the curve is proportional to the fre- quency of the waves. If we increase the frequency of the waves (radiation) the radius gets larger and the curve straightens up. It can be said that the slice theorem is a good approxima- tion of the diffraction theorem for high-frequency radiation sources like X-rays [10]. The diffraction theorem is used in a similar way in the reconstruction of the images as the slice theorem in the backprojection method in the frequency domain (compare to the slice theorem in eq. (1.5)). 1.3.3 Reflectivity imaging The pioneering work on this topic was done by Norton and Linzer [52]. "Reflectivity" is understood as a property of the tissue causing the change of direction of the propagating ultrasound signal. It is not only the scattering but also refraction on the boundaries of parts of tissue with different acoustic impedance. The theory was originally designed for a circular two-dimensional scanning geometry, but it is simply extendable to three dimensions. The previously mentioned methods make it possible to reconstruct the local attenuation or speed of sound properties of the scanned object. In velocity imaging it is the time-of-flight, in attenuation imaging it is the shape and magnitude of the first pulse in the recorded signal, which are used as the source of data for the reconstruction. Both modalities only utilize the properties of the first pulse giving information about the object along the straight line (the shortest path) between the emitter and receiver. Reflectivity imaging, however, exploits the whole length of the recorded signals to ob- tain information about the object from non straight paths of propagation of the ultrasonic waves. The method somewhat resembles the B-mode ultrasound in that it uses the backscat- tered waves as the input data to build the image. It also goes further, not only because it is able to combine B-mode images from all sides of the object, but also use information of the scatter- ing to all sides (not just the backscattering). The reconstruction algorithm builds the image in the following way: each point in the resulting image accumulates values from the recorded signals, corresponding to the length of the reflection paths, that is, the sum of the paths from the emitter to the reflection point and from the reflection point to the receiver. Figure 1-4: Reconstruction principle of reflectivity based on the assumption that the sound speed is constant (or similar enough) inside and outside of the scanned object. (Taken from [65]) The method assumes that the speed of sound is constant in and outside of the scanned object (Figure 1-4). A modification of the method assumes first reconstructing the velocity 12 map of the object and then using this information for the reflectivity imaging to account for the different speeds of sound. Although no physical properties are reconstructed in reflectivity imaging, the images usually have a high information value considering the structure of the imaged object. Moreo- ver, the images can be reconstructed with a sub-millimeter resolution. 1.3.4 Algebraic Reconstruction Techniques Besides the backprojection algorithm there is another category of methods solving the problem of reconstructing images from projections. The algebraic methods are based on dis- cretization of the projection area. The continuous function 𝑓(𝑥,𝑦) is approximated by means of its samples 𝑓𝑖 ,𝑘 via an interpolation 𝑓 𝑥,𝑦 = 𝑓𝑖 ,𝑘 𝑔𝑖,𝑘 𝐾 𝑘=1 (𝑥, 𝑦) 𝐼 𝑖=1 (1.7) The functions 𝑔𝑖 ,𝑘 are usually staircase functions (having a constant non-zero value in the area of the pixel 𝑖,𝑘), but generally can be any higher order interpolation function. It is possible to rearrange 𝑓𝑖 ,𝑘 into a column vector f and create a matrix W of weights with each row corresponding to the intersections of one beam with the scene pixels. Then we can express the set of measured values as the vector 𝐖𝐟 = 𝐩 (1.8) This approach has the advantage that the ray-integrals in the projections do not neces- sarily need to be parallel or even straight! The shape of each ray (beam) is encoded into each row of the weight matrix W. To solve for the original image, the system only needs to be inverted 𝐟 = 𝐖−1𝐩 (1.9) And so, the problem of reconstructing an image has been transformed into a problem of solv- ing a system of linear equations. Note that the matrix W does not need o be square (depending on the number of projections) and in that case the resulting image is a least mean square solu- tion (by e.g. pseudoinversion) of the problem. Because the system of equations is usually too large to be solved by the conventional approaches (the number of elements of W can reach up to the order 10 8 to 10 12 ), an iterative solution (SART, SIRT, Newton-Kazmarz method [28],[29]) is used. 1.4 The need for calibration For the reconstruction of the tomographic images it is crucial to know the properties of the imaging system and especially the used transducers. For example in attenuation imaging, the frequency content of the first pulse must be correctly determined in order to truthfully re- construct the local distribution of attenuation in the object. The system usually operates at frequencies near the transducers‟ resonance and the transfer function of the transducers is usually highly variable up or down the spectrum. The frequently used wideband ultrasonic pulses get distorted first in the emitter and then also in the receiver transducer. After such distortions, the exact frequency content of the pulse is hard to determine correctly and some calibration is necessary. The transducers also usually have a strong angular dependence. Besides a wide main lobe, side lobes can also be present which (if not accounted for in a calibration step) can cause strong artifacts in the reconstructed images. Both of the above described phenomena (the frequency and angular dependencies of the transducers) are usually known already during design time, and can therefore be built into some internal correction of the system from the beginning. But transducers are subject to aging and these properties change. Especially the overall sensitivity might degrade by several magni- 13 tudes. Such transducers add a lot of noise to the signals and the reconstruction, without a prop- er calibration, is again prone to give faulty images full of artifacts. Finally, new USCT systems consist of thousands of independent transducers and it may not always be possible to build them with their positioning so precise as is needed for the re- construction of e.g. velocity images. Differences of tenths of millimeters already cause devia- tions of several meters per second when estimating the speed of sound. In reflectivity images, small positioning errors might cause phase cancelation and again obscure the reconstructed image in artifacts. The calibration of USCT systems is therefore a very important aspect of image recon- struction which, if applied properly, can lead to significant image quality improvement. 14 2 State of the art in USCT This chapter discusses the recent advances in science in the field of ultrasound comput- er tomography and calibration techniques. First, a list of present-day attempts to develop a fully functional ultrasound computed tomography system is presented in chapter 2.1. Only a brief overview is devoted to each project. More space is dedicated to the partnering Karlsruhe project (chapter 2.2), in which the author took part. Chapters 2.3 and 2.4 are over viewing the up-to-date methods and techniques of trans- ducer sensitivity calibration and position calibration. The presented methods are used in cur- rent ultrasound computer tomography systems, but also in other technical areas and were taken as a basis in development of the new calibration techniques by the author. 2.1 Ultrasound computed tomography systems 2.1.1 High-Resolution Ultrasonic Transmission Tomography System A very promising new ultrasound system was described in IEEE Transactions on Med- ical Imaging in 2005. This “High-Resolution Ultrasonic Transmission Tomography System” (HUTT) was developed in the Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, USA by Jeong et.al. [30]. Contrary to most of the other ultrasound tomography systems, this system uses high- frequency transmitters (center frequency 8 MHz with 50% bandwidth). The transducer ele- ments are only 0.4mm x 0.4mm. The received signals are sampled at 100 MHz in a 14 bit A/D converter. The part of the recorded signal which contains the first-arrival pulse is extracted and processed for multi-band analysis, utilizing the frequency dependent characteristics of the acoustic attenuation. The images are built by fusing conventional backprojection sinograms, each of which is set up for one of the frequency bands. The fusion is realized by a Local Principal Compo- nent Analysis that maps all of the narrow-band 3D sinograms into one, which is then recon- structed using the classical backprojection technique and a Shepp-Logan filter. Figure 2-1: Reconstruction of a sheep kidney. The images were coregistered from three modalities: optical (left), MRI (middle) and HUTT (right). (Taken from [30]) 2.1.2 The Ring Transducer System for Medical Ultrasound Research This 2D ultrasound computed tomography system (Figure 2-2) has been developed for experimental studies of scattering and imaging. The ring transducer array (built by Nihon Dempa in Kogyo Co., ltd., Tokyo, Japan) consists of 2048 rectangular elements with a 2.5- MHz center frequency, a 67% - 6 dB bandwidth, and a 0.23-mm pitch arranged in a 150-mm- diameter ring with a 25-mm elevation. At the center frequency, the element size is 0.30  x 42  and the pitch is 0.38 [73]. 15 Figure 2-2: Ring transducer system photo. The ring transducer is in the center of the cylindrical apparatus and forms a portion of a cylindrical water tank that is temperature-controlled to 30◦C ± 0.2◦C. The electronics are in the background. A motorized gantry located above the ring permits test objects to be raised or lowered through the plane of the ring. (Taken from [73]) The system has 128 parallel transmit channels, 16 parallel receive channels, a 2048:128 transmit multiplexer, a 2048:16 receive multiplexer, independently programmable transmit waveforms with 8-bit resolution, and receive amplifiers with time variable gain independently programmable over a 40-dB range. Receive signals are sampled at 20 MHz with 12-bit resolu- tion. Arbitrary transmit and receive apertures can be synthesized. This system also incorporates a calibration mechanism to compensate for unavoidable element-to element variations in sensitivity and time response and from deviation in element position from an ideal circle. The used algorithm is described in detail in chapter 2.4.3. 2.1.3 A diffraction tomography system This system has been developed by Andre et. al. (Department of Radiology, Veterans Affairs Medical Center, San Diego, CA) for the purposes of experimenting with ultrasound computed tomography and its applications to breast imaging. Low-power discrete frequency sound in the range of 0.3–1.2 MHz, two cylindrical arrays of 512 and 1024 PZT transducers, and high spatial sampling of the wavefront are used [4]. As for the reconstruction, the system uses a diffraction tomographic reconstruction me- thod. One transducer at a time is activated and allowed to reach steady state at which point the remaining transducers measure the phase and amplitude of the ultrasound signal. The image is formed by inverting the wave equation and calculating a complex scattering potential at all locations in the 2D plane. The wave equation simplifies to an expression with Henkel and Bessel functions because the transducers work only in discrete frequencies. 16 Figure 2-3: Nine sequential image slices through entire breast with ~4-mm overlap from nipple to the chest wall. Magnitude with aberration correction. (Taken from [4]) The team published very promising reconstruction images of in-vivo breast as early 1997 [2][3][4], and for a few following years had no activity. More recently an analysis of patient breast images from a clinical trials series of 25 patients has been published [33]. The study suggests that ultrasound diffraction tomography has the potential for discriminating be- tween different tissue types. The authors also conclude that the image contrast values of ma- lignant lesions provide sufficient discrimination from normal breast tissue to be potentially of diagnostic value. 2.1.4 Computerized Ultrasound Risk Evaluation (CURE) program The Karmanos Cancer Institute (KCI) in collaboration with Lawrence Livermore Na- tional Laboratory (LLNL) has built an ultrasound scanning system designed to simulate vir- tually any transducer array design [40]. The system has six degrees of freedom and has been optimized for signal-to-noise (SNR), pulse shape, minimization of systematic errors, and au- tomated preprocessing. As the authors claim, the system is very flexible. In one of the documented experi- ments, the system comprised of two identical line transducers (0.38 mm wide and 12 mm high) movable in a circular geometry. The relative positions of the two transducers positioned with accuracy 0.05 mm or higher. For each position of the transducer, the receiver was moved around a 320 degree circular path of radius 15 cm. The transmitter emitted a pulse for each position of the receiver for up to 1600 receiver positions (position increments of 0.2 deg). The transmitter position was also moved along the circular path and the firing sequence repeated for each new position of the transmitter [16]. The targets were placed at the center of the ring with the long axis of the cylindrical phantom oriented vertically relative to the plane of the ring. Each data set represents a 2-D slice through the target. The ring plane was translated in the vertical direction allowing for 3-D reconstructions from stacked 2-D planes of data. All scans were performed at 10 mm slice thickness, as determined by the beam width of the transducers. The data sets resulting from the above scans were typically ~ 2 GB per slice. The raw data sets consisted of measured pulse trains, processed to determine the frequency spectrum and arrival time for each pulse. The arrival times were used to construct reflectivity and sound speed images while the pulse spec- trum data were used to determine the attenuation [16]. 17 The collaborating groups tried several different reconstruction approaches resulting in a comparison study published in [40]. An example of a reconstructed image can be seen in Figure 2-4. Figure 2-4: Full aperture tomography (FAT) reconstruction of a cadaveric breast placed in formalin and sealed in a 100mm container. The phantom surrounded by the transducer arms of the scanner can be seen on the right. (Taken from [40]) More recently, the group has carried out a study with 19 patients whose mammograms and follow-up ultrasound identified suspicious masses. The CURE exam was interposed be- tween the standard US exam and subsequent biopsy. Biopsy results were therefore available for all 19 patients. Based on the preliminary CURE data the group has utilized six CURE di- agnostic criteria for cancer. In the small sample, it appears that women with higher scores are more likely to have cancerous masses [17]. 2.1.5 Karlsruhe USCT system Finally, the USCT system built in the Forschungszentrum Karlsruhe (FZK) will be giv- en a separate chapter as all of the data which the author worked on originated from here. 2.2 Karlsruhe USCT system and previous results The initiator of the whole project, Dr. Stotzka, developed the first experimental 2D prototype of the USCT system in 2001 [65]. Since then, the topic has been extensively studied in FZK [50][57][66][67] and today a new experimental 3D system is functional [68]. Because of a conceptual difference in the construction of the 2D and the 3D systems a chapter will be devoted to each separately. 2.2.1 The Karlsruhe 2D USCT system The 2D system developed in FZK was designed to cover a ring around the imaged ob- ject by as many transducers as possible. To keep the system reasonably simple, two moving arms are used. Each of the arms carries an ultrasonic probe with 16 transducers and could be circled around a steel ring with 100 fixed positions. In the first probe group only one of the 16 transducer elements is used for transmitting ultrasonic pulse waves. All 16 elements of the second probe group are used for receiving. The individually received signal, so called A-scan, 18 is recorded by the receiving element each time the sending element is excited by an electrical pulse. The whole system scan is made in 100 steps. Each step consists in successively firing the emitter element at one position and changing the receiving probe position around the ring while recording the received A-scans. In the next step, the sending element is moved into the next position and the receiving probe is again circled around the ring. At each step 1600 A-scans are recorded, therefore a full system scan consists in 160,000 A-scans. The ring has a diameter of 12 cm. Each transducer element is 0.2 mm wide and 10 mm high. The array of the 16 elements has a pitch of 0.25 mm. The center frequency of the ultra- sonic pulses 2.5 MHz was chosen as a compromise between large absorption at higher fre- quencies and lower resolution due to larger wavelength at lower frequencies. To achieve very short pulses, the transducers are strongly damped and transmit broadband pulses. The received 200 s long signals are digitized with a sampling rate of 50 MHz and signal quantization of 10 bits. Figure 2-5: Experimental prototype of the 2D USCT system in Forschungszentrum Karlsruhe. The white tank is filled with water as a coupling medium for the examined object. Two ultrasonic probes are mounted on ro- tated rings to simulate all emitter and receiver positions. (Taken from [65]) 2.2.2 The Karlsruhe 3D USCT system The 3D USCT system, besides adding a new dimension, was designed to overcome some of the difficulties of the 2D system. Mainly the data acquisition time was overwhelming on the 2D system, because the moving probe arms had to be moved manually and each A-scan was individually recorded by a digital oscilloscope. The 3D system on the other hand handles the acquisition automatically with a PC. There are also many more physical transducer ele- ments, which can simultaneously record A-scans. 2.2.2.1 System architecture The new 3D ultrasound computer tomography system consists of three parts: a tank containing the sensor system, data acquisition hardware and a computer workstation (Figure 2-6). The tank has a diameter of 18 cm and a height of 15 cm. 48 transducer array systems (TAS) are mounted into the tank walls carrying each 32 receiving and 8 emitting elements (Figure 2-7). The transducer elements can be accessed individually. The resulting cylindrical array is rotated in 6 steps to fill the gaps between the transducers. Thus a fully covered cylin- drical array can be emulated, resulting in a total of 9216 receiver and 2304 emitter positions [68]. Because of the rotation, the A-scans can be sent from and received at virtually any po- sition on the cylinder, but not all emitter-receiver position combinations are possible. Still the 19 total number of A-scans produced by a full system scan is approximately 3.5 million. Al- though this number is higher than the total number of A-scans in the 2D system, there are less A-scans per a z-layer in the 3D system (as there is only one layer in the 2D system). The received signals and the control signals are gathered in four blocks containing preamplifiers, address generation and control logic. Figure 2-6: Overview of the 3D USCT system in Forschungszentrum Karlsruhe. Left: the cylindrical tank with 48 mounted transducer array systems. Middle: Data acquisition unit with 192 channels. Right: computer sys- tem for image reconstruction. (Taken from [69]). The data acquisition electronics is a modified design of the system for the Auger fluo- rescence detector [19]. It is based on 9 9HE-card boards connected by a modified VME-bus. One board controls the data acquisition process, the data storage and the transfer to the com- puter workstation. The other 8 boards are carrying 24 data recording channels each. A channel consists of the analog signal processing, an A/D converter (10 MHz, 12 Bit) and the digital signal processing. The digital processing is based on an array of four FPGAs executing the online data storage, noise reduction and simple data reduction [68]. While one emitting transducer is activated with coded excitation, all receiving trans- ducers receive simultaneously the scattered signals. 192 channels are sampled and recorded in parallel. 48 multiplexing steps are needed to record an A-scan from every receiver. The data is transferred to a computer workstation and the next transducer-emitter is selected. The data storage, the image reconstruction and the visualization of the results are accomplished on the computer workstation [68]. The recorded A-scan signals are sampled at 10MHz and 12bits. A complete scan pro- duces about 20GB of data. Although theoretically it is possible to scan the whole system with- in a few seconds, it presently takes 10+ hours. This is because each A-scan is a result of aver- aging about 10 to 50 measurements to gain a higher signal to noise ratio. The measurement also has to be multiplexed in time (only 192 of 1536 receiver element signals are recorded at a time – using all available channels) to allow for the relatively slow process of storing the data on a PC hard drive. In the near future, it is planned to install a new data acquisition unit, which will allow for scan times under 10 minutes. 2.2.2.2 Transducer array systems - TAS The ultrasonic transducers are grouped into blocks called transducer array systems (TAS). The transducer elements have the mean frequency of about 3 MHz. They are made of a piezzo-ceramics plate, which is structured by a sub-dicing technique [64], to produce an ele- ment of 1.4 mm by 1.4 mm. This results in a relatively narrow angular characteristic. The de- sign resulted from a compromise between an ideally omnidirectional characteristic and the need to focus the ultrasonic energy into the center of the USCT tank to get the waves through the examined object. 20 Figure 2-7: The transducer array system (TAS) – folded and sealed (left) and an interior view of an unfolded TAS (right). The transducer array systems are designed and manufactured in a hybrid laboratory of the FZK. The almost entirely automatic manufacture process guarantees high quality and high reproducibility of the transducer characteristics (±2%) and low costs (~100 € per transducer array system including front-end electronics) [68]. The ultrasonic pulses generated by the TASes are controlled by the so called coded ex- citation pulse – a short electrical signal generated by a D/A converter on the data acquisition boards and programmable by the controlling software. The ultrasonic pulses are therefore con- trollable in terms of length, frequency, bandwidth, and so on. The power spectrum of the transmitted ultrasonic pulses is of course highly dependent on the transfer functions of the electronics and the ultrasonic transducers. 2.2.3 FZK Reconstruction method (the sum-and-delay algorithm) The following three imaging modalities are possible using the Karlsruhe USCT: reflec- tivity imaging, speed of sound imaging, and attenuation imaging. Since the introduction of the 3D system (2004), the group mainly concentrates on the reflectivity imaging, because the transducer array coverage is very sparse and achieves full coverage only by rotating the whole cylinder with both emitter and receiver elements. The algorithm used by the Karlsruhe team for the reflectivity imaging is a so-called sum-and-delay algorithm [70]. For each point in the image, the amplitude (or some prepro- cessed version) of the acquired A–scans at the position corresponding to the distance between emitter and receiver, is accumulated 𝑓 𝑥 = 𝑇 𝐴 𝑗 ,𝑘 𝑎𝑗 + 𝑏𝑘 𝑐 𝑗 ,𝑘 (2.1) where f denotes the reflectivity image, 𝑥 the coordinates of the reconstructed point, T are pre- processing steps, A(j,k) the A–scan acquired at sending position 𝑥 𝑗 and receiving position 𝑥 𝑘 . c is the speed of sound in water, and aj and bk are the distances from the reconstructed point to the emitter and receiver, respectively (also see Figure 1-4). Preprocessing steps T can be also implemented e.g. matched filtering, envelop calculation, deconvolution, etc. Equation (2.1) is valid for constant speed of sound, small attenuation, spherical emittance and Huygens‟s point scatterers [58]. Using the same assumptions and neglecting noise, the resulting value of each voxel can be given as: 𝑓 𝑥 = 𝑅 𝑙 𝑑𝑙 = 𝑁𝑅 𝐸 𝑖 ,𝑘 𝑥 𝑗 ,𝑘 + 𝜀 𝐼,𝐾 𝑥 (2.2) where is 𝑅 𝑥 the reflectivity map of the imaged object, E(i,k) is the integral over the surface of the ellipsoid given by the emitter i, the receiver k, and voxel position 𝑥 𝑘 . Then the resulting voxel amplitude 𝑓 𝑥 is the sum of the reflectivity R at 𝑥 magnified by N, the number of ap- plied emitter–receiver combinations, and an error term 𝜀(𝐼,𝐾) 𝑥 (I, K indicate the 21 dependence of ε from the geometry and discretization of the aperture). ε is the result of mul- tiple scatterers received at the same time instance [58]. To test the imaging capabilities of the 2D system, the team in Karlsruhe built a phan- tom consisting of several tubes with gelantine, straws, and threads. The smallest structures within the phantom are nylon threads with a diameter of 0.1 mm each, corresponding to ap- proximately 0.2 wavelengths of the emitted ultrasound signal [68]. Figure 2-8: Phantom and reconstructed images from Karlsruhe 2D USCT. Left: rough blueprint of the phan- tom. The smallest structures consist of nylon threads with a diameter of 0.1 mm and a spacing of 0.5 mm. Right: reconstructed cross-section. (Taken from [68]) To mimic the imaging of a real breast, a clinical breast phantom (triple biopsy breast phantom, CIRS, Inc., Norfolk, USA) was imaged with the Karlsruhe 3D USCT system. The phantom contained cancer and cyst mimicking masses of 2 to 10 mm in diameter. For compar- ison a MR volume of the phantom (1.5T Siemens Magnetom Vision, double breast coil, T1– weighted, (1.37 mm 3 voxel size) was acquired (Figure 2-9, left). Figure 2-9:.Comparison of MRI and Karlsruhe 3D USCT images of a breast phantom. Left: an MR image. Middle: USCT image with artifact reflections from a “cancerous” region with high reflectivity. Right: cor- rected USCT image with the artifact removed and normalized reflectivity values (the smaller magnifications show the two partly obscured “cysts”). (Taken from [58]) The 3D USCT reflectivity image reconstructions are visible in (Figure 2-9, middle and right). As ultrasound images the tissue borders, the USCT images are expected to be similar to the gradient images of the MRI. The ”skin” of the phantom is clearly visible, also one “cyst” (round dark area) and the “cancer” (partly white double circle) structure at the top left is im- aged with well defined outlines. A bright blob approximately in the center of the image ob- scures the other two “cysts”. This artifact is caused by a cluster of strong reflectors (”cancer”) 22 near the “nipple” of the phantom outside the focal region. For an image based on A–scans where the reflections from this point have been deleted; the bright artifact vanishes and the obscured “cysts” become partly visible. 2.3 Published sensitivity calibration methods The transducer array elements used in any ultrasonic system can never be manufac- tured with 100% repeatability and differences in quality are unavoidable. Each transducer can be described with its own radiation patter – an angle and frequency dependent function de- scribing the emitted pressure field. The transducer to transducer variations in the shape of this function are usually negligible, but the differences in the magnitude of this function can be considerable. If these variations are compensated for, the quality of the reconstructed images can be significantly improved. Systems with only a small number of ultrasonic transducers can utilize the classical ap- proach of measuring the radiation field using a calibrated hydrophone. The measurement can be done for each element individually yielding a sensitivity parameter or even the whole radia- tion pattern. The received signals are then multiplied by the corresponding correction factors. Such an approach was used by Andre et.al. in their diffraction tomography system [4]. A different approach to sensitivity calibration was used in the Ring Transducer System for Medical Ultrasound Research (overview of the system can be found in chapter 2.1). The calibration takes advantage of the fact that the individual transducers act as both the emitters and receivers. The calibration uses a 0.10-mm-diameter metal wire at the approximate center of the ring shaped transducer array, as depicted in Figure 2-10. A single element is used to transmit a pulse and then receive the echo from the wire. This pulse-echo experiment is repeated 10 times for each of the 2048 elements and the repetitions at each element are averaged [73]. Figure 2-10: Calibration geometry. For calibration, a pulse is transmitted from a single element, and the echo from a wire positioned at the center of the ring is received by the same element. (Taken from [73]) A reference waveform 𝑠 (𝑡), found by averaging all of the time-aligned waveforms, is formed first. The reference waveform is then used with the pulse-echo signal 𝑠𝑖(𝑡), where i is the element number, to compute finite-impulse response compensation filters 𝑟𝑖(𝑡) that match the pulse-echo signal from the wire with the reference waveform. The computation is based on the relation: 𝑠𝑖(𝑡) = 𝑒 𝑡 ⊗ 𝑕𝑖(𝑡) ⊗ 𝑕𝑖(𝑡) (2.3) where 𝑕𝑖(𝑡) is the impulse response of element i, 𝑒 𝑡 is the excitation signal applied to the element, and ⊗ represents temporal convolution. The relation for the calibration filter 𝑟𝑖(𝑡) applied before transmission and after reception is then defined by: 23 𝑠 (𝑡) = 𝑒 𝑡 ⊗ 𝑟𝑖(𝑡) ⊗𝑕𝑖(𝑡) ⊗ 𝑕𝑖(𝑡) ⊗ 𝑟𝑖(𝑡) (2.4) in which the pulse-echo response at element i has been set equal to the average signal. Substi- tuting (3) in (4) and solving for the calibration filter at each element yields: 𝑅𝑖 𝑓 = 𝑆 (𝑓) 𝑆𝑖(𝑓) (2.5) where all the uppercase symbols denote the temporal Fourier transform. The values of 𝑅𝑖 in the system passband are the coefficients of the corresponding finite-impulse-response com- pensation filter [73]. Figure 2-11: Average RF echo waveform (image on top) from a wire target and the corresponding power spec- trum (center) and phase (bottom). These are the reference waveform and spectrum used for the computation of the calibration filters. (Taken from [73]) Such an approach elegantly solves not only the scalar sensitivities of individual trans- ducer elements but also gives a temporal compensation filter. On the other hand, the directivi- ty pattern is omitted completely in the calculations and so such a filter is only valid for omni- directional transducers. The above described calibration approach shows how an ultrasound system with thou- sands of transducers can be automatically calibrated without using any additional measuring equipment. Unfortunately, this particular method could not be used in case of the Karlsruhe USCT mainly because of the lack of transducer elements, which are able to both emit and re- ceive. 2.4 Known position calibration methods The following chapters give an overview of several position calibration techniques, which are used in other technical areas and applications. It is shown and argued why none of these methods are fully suitable for the case of the Karlsruhe USCT position calibration and that therefore a novel method had to be developed. 2.4.1 Global Positioning System (GPS) The Global Positioning System (GPS), although primarily a military application devel- oped by the US Army, has been promoted for civil uses since the beginning of its develop- ment, and has therefore been very well described in numerous scientific and technical publica- tions. The GPS was conceived as a ranging system from known positions of satellites in space to unknown positions on land, at sea, in air and space. The original objectives of the system 24 were the instantaneous determination of position and velocity and the precise coordination of time. The system consists of three major segments [26]:  the space segment consisting of satellites orbiting the Earth which broadcast the position- ing signals,  the control segment steering the whole system,  the user segment including many types of receivers. Because many aspects, which are dealt with in GPS (satellite orbits, signal encoding, atmospheric effects, relativistic effects, etc.), are not necessary to deal with in the USCT, only the basic mathematical concepts will be discussed here. Despite the large number of literature sources on GPS the following paragraphs dis- cussing the main concepts are mostly taken from [26] as the author found it as the most com- prehensive. 2.4.1.1 Pseudorange equations The GPS uses pseudoranges (pseudo distances) resulting from receiving and analyzing the broadcast satellite signal. The pseudorange is derived from measuring the travel time of the coded signal from the satellite to the receiver t and multiplying it by its velocity c (the speed of light). The clocks of the satellite and the receiver are never perfectly synchronized. Instead of true ranges , pseudoranges R are obtained where the synchronization error  (denoted as clock error) is taken into account (2.6). Ri j = ρi j + cΔδi j (2.6) where the lower indices refer to the receiver site and the upper indices refer to the satellite number in view of the receiver. Let us now take a look at the unknowns in equation (2.6). The distance from the satel- lite to the receiver can be explicitly written as: ρi j = (Xj − Xi)2 + (Yj − Yi)2 + (Zj − Zi)2 (2.7) where {X j , Y j , Z j } is the known position of the satellite (coded in the broadcast signal) and {X- i, Yi, Zi} is the unknown position of the receiver in given absolute (world) coordinates. The clock error i j is composed of two parts: the unknown bias i of the receiver clock from the universal GPS time and the known bias j of the satellite from the GPS time, which is also encoded in to the broadcast signal. Δδi j = δi − δ j (2.8) Consequently each equation (2.6) comprises of four unknowns: the three point coordi- nates of the receiver unit, and the clock bias. Thus, four of these equations (and therefore four satellites in view of the receiver unit) are always necessary to solve for the four unknowns. Substituting (2.8) into (2.6) and shifting the satellite clock bias to the left side of the equation yields: Ri j + cδ j = ρi j + cδi , ∀j (2.9) Note that the left side of the equation contains observed or known quantities, while the terms on the right side are unknown. 2.4.1.2 Linearization Note that the true distance  in (2.7) is nonlinear with respect to the unknown position coordinates of the receiver. In order to solve for the unknowns a linearization step is useful. Let us assume approximate values Xi0, Yi0, Zi0 for the unknowns, an approximate dis- tance ji0 can be calculated by 25 ρi0 j = (Xj − Xi0)2 + (Yj − Yi0)2 + (Zj − Zi0)2 ≡ f(Xi0, Yi0, Zi0) (2.10) Using the approximate values, the unknowns Xi, Yi, Zi can be decomposed as Xi = Xi0 + ΔXi Yi = Yi0 + ΔYi Zi = Zi0 + ΔZi (2.11) where now Xi, Yi, Zi are the new unknowns. This means that the original unknowns have been split into a known part (represented by the approximate values Xi0, Yi0, Zi0) and an un- known part (represented by Xi, Yi, Zi). The advantage of this splitting-up is that the func- tion f(Xi, Yi, Zi) is replaced by an equivalent function f(Xi0+Xi, Yi0+Yi, Zi0+Zi) which can be now expanded into a Taylor series around the approximate point. This leads to f(Xi, Yi, Zi) ≡ f(Xi0 + ΔXi, Yi0 + ΔYi, Zi0 + ΔZi) = f(Xi0, Yi0, Zi0) + ∂f(Xi0, Yi0, Zi0) ∂Xi0 ΔXi + ∂f(Xi0, Yi0, Zi0) ∂Yi0 ΔYi + ∂f(Xi0, Yi0, Zi0) ∂Zi0 ΔZi+. . . (2.12) where the expansion can be truncated after the linear term thanks to the closeness of the esti- mate point to the exact position; otherwise the unknowns would appear in nonlinear form again. The partial derivatives are obtained from (2.10) by ∂f(Xi0, Yi0, Zi0) ∂Xi0 = − Xj − Xi0 ρ i0 j ∂f(Xi0, Yi0, Zi0) ∂Yi0 = − Yj − Yi0 ρ i0 j ∂f(Xi0, Yi0, Zi0) ∂Zi0 = − Zj − Zi0 ρ i0 j (2.13) and are the components of the unit vector pointing from the satellite towards the approximate site. The substitution of equations (2.10) and (2.13) into (2.12) gives ρi j = ρi0 j − Xj − Xi0 ρ i0 j ΔXi − Yj − Yi0 ρ i0 j ΔYi − Zj − Zi0 ρ i0 j ΔZi (2.14) where the equivalence of f(Xi, Yi, Zi) with ji has been used. This equation is now linear with respect to the unknowns Xi, Yi, Zi. 2.4.1.3 Solving the system of equations Substituting equation (2.12) into (2.9) yields a linearized pseudorange equation with known values on the left and unknown on the right Ri j − ρi0 j + cδ j = − Xj − Xi0 ρ i0 j ΔXi − Yj − Yi0 ρ i0 j ΔYi − Zj − Zi0 ρ i0 j ΔZi + cδi (2.15) Because four variables are unknown, four equations have to be set up to have the solu- tion determined. The shorthand notations 26 l j = Ri j − ρi0 j + cδ j aXi j = − Xj − Xi0 ρ i0 j aYi j = − Yj − Yi0 ρ i0 j aZi j = − Zj − Zi0 ρ i0 j (2.16) help to simplify the representation of the system of equations. Assuming now four satellites numbered from 1 to 4, then l 1 = aXi 1 ΔXi + aYi 1 ΔYi + aZi 1 ΔZi + cδi l 2 = aXi 2 ΔXi + aYi 2 ΔYi + aZi 2 ΔZi + cδi l 3 = aXi 3 ΔXi + aYi 3 ΔYi + aZi 3 ΔZi + cδi l 4 = aXi 4 ΔXi + aYi 4 ΔYi + aZi 4 ΔZi + cδi (2.17) is the appropriate system of equations. Note that the superscripts are the satellite numbers and not exponents! Introducing A = aXi 1 aYi 1 aZi 1 c aXi 2 aYi 2 aZi 2 c aXi 3 aYi 3 aZi 3 c aXi 4 aYi 4 aZi 4 c x = ΔXi ΔYi ΔZi δi l = l 1 l 2 l 3 l 4 (2.18) the set of linear equations can be written in the matrix-vector notation l = Ax (2.19) For this first example of a linearized GPS model, the re-substitution of the vector l and the matrix A using (2.16) is given explicitly by: l = Ri 1 − ρi0 1 + cδ1 Ri 2 − ρi0 2 + cδ2 Ri 3 − ρi0 3 + cδ3 Ri 4 − ρi0 4 + cδ4 A = − X1 − Xi0 ρi0 1 − Y1 − Yi0 ρi0 1 − Z1 − Zi0 ρi0 1 c − X2 − Xi0 ρi0 2 − Y2 − Yi0 ρi0 2 − Z2 − Zi0 ρi0 2 c − X3 − Xi0 ρi0 3 − Y3 − Yi0 ρi0 3 − Z3 − Zi0 ρi0 3 c − X4 − Xi0 ρi0 4 − Y4 − Yi0 ρi0 4 − Z4 − Zi0 ρi0 4 c (2.20) The coordinate differences Xi, Yi, Zi and the receiver clock error i result from the linear system (2.19) x = (A T A)−1A T l (2.21) (here the pseudoinverse LMS solution is used, so that the matrix A doesn‟t have to be square as in the case, when more than four satellites are observable from the receiver site) 27 The desired point coordinates Xi, Yi, Zi are then finally obtained by updating the coor- dinate estimates with the calculated error values (2.11). Although it is not usually stressed enough in the GPS literature, one or more linearize- solve-update iteration steps may be needed to reach the desired position accuracy. Note that the above iterative approach is actually an implementation of the Gauss-Newton method [56] to solve the original nonlinear equation system. First, a vector of estimate values Xi0, Yi0, Zi0 is set up. Then, from the nonlinear set of equations (2.9), a Jacobian matrix A is computed. The next step comprises of solving the li- nearized equation set (2.19). And finally the vector of estimates is updated by the solved vec- tor of error values (2.11). These steps can be iterated until the root-mean-square of the solved error vector falls below a predetermined threshold. 2.4.1.4 Dilution of Precision (DOP) An important factor in achieving high quality positioning results with GPS (and any other positioning system) is a low dilution of precision (DOP). The DOP is a measure of the relative geometry of the visible satellites and greatly affects the resulting position uncertainty. A simple example shows the effect of the geometry of two transmitters in a basic 2D positioning system (Figure 2-12). In a) the transmitters form nearly a right angle with respect to the receiver, and the uncertainty of the distance measurements (indicated by the patterned areas) results in a small position uncertainty (low dilution-of-precision). In b) the transmitters form a sharp angle and the position uncertainty of the receiver is considerably larger (high dilution-of-precision). Figure 2-12: DOP in a basic positioning system (taken from [37]). If we look at how the receiver position is obtained from the measured pseudoranges, we can see that the vector of estimate-error terms x in (2.19) has to be solved by pseudo- inverting the system matrix A x = (A T A)−1A T l (2.22) It is assumed that all pseudorange measurements are equally uncertain and that no correlations exist between the errors. Otherwise a more general form can be used: x = (A T WA)−1A T Wl (2.23) where W is a weight matrix reflecting the differences in errors of the measurements and the correlations among them. This weight matrix is also equal to ς0 2CΔPC −1 , in which CΔPC is a covariance matrix of the pseudorange errors and ς0 2 is known as the prior variance of the unit length. In general, the solution of a nonlinear problem must be achieved by iteration to obtain the result. However, if 28 the linearization point is sufficiently close to the true solution, then only one iteration is re- quired [37]. If we are asking how accurate can the position and time-delay parameters be calculated for a GPS receiver, we are actually asking how do the pseudorange measurements and model errors affect the estimated parameters in (2.23)? This is given by the covariance law: CΔx = [(A T WA)−1A T W]CΔPC [(A T WA)−1A T W]T = (A T CΔPC −1 A)−1 (2.24) in which CΔxis the covariance matrix of the parameter estimates. If we now assume that the measurement and model errors are Gaussian and with the same standard deviation  for all observations, and that they are uncorrelated, then CΔPC = Iς2, where I is the identity matrix. In this case equation (2.24) simplifies to: CΔx = (A T A)−1ς2 = Dς2 (2.25) With a value for , we can compute the components of Cx using the above equation. We then can get the measure of the overall quality of the least-squares solution by taking the square root of the sum of the parameter estimate variances: ςG = ςE 2 + ςN 2 + ςU 2 + ςT 2 = ς D11 + D22 + D33 + D44 (2.26) in which E  , N  , and U  are the variances of east, north, and up components of the receiver position estimate, and T  is the variance of the receiver clock offset estimate. The above is an estimate of the solution accuracy. It is dependent on the measurement and model error standard deviation () on one hand, and the square-root of the trace of matrix D on the other. The scaling factor formed from the components of matrix D is only dependent on the satellite geometry. It is usually greater than one and thus it amplifies the pseudorange error, or dilutes the precision, of the position determination. This scaling factor is therefore usually called the geometrical dilution of precision: DOP = trace(D) (2.27) 2.4.1.5 The use of GPS principles in USCT This chapter covered the principles of the GPS positioning method. As it is, the GPS method cannot be used for anything else than for what it was designed: navigation. The impor- tance for the use in USCT lies in the fact that it is possible to obtain positions of the receiver only by observing the time-of-flight of the signals travelling from the orbiting satellites to the receiver. Also the mathematical apparatus of the dilution of precision proved useful in the novel USCT positioning method, where the author used it to optimize the anchoring (boundary conditions) of the system of equations. 2.4.2 A published GPS modification for an ultrasonic systems In [41], Yue Li proposed a novel calibration method based on the principles of the GPS positioning. The author needed to calibrate the positions of ultrasonic transducers, which formed a large underwater imaging system. The ultrasonic receivers substituted GPS receivers, and a high-precision positioning device with a hydrophone replaced the GPS satellites. The method was also extended to include the calibration of the speed of sound of the water in which the imaging system was submerged. The imaging system comprises one powerful ultrasonic transmitter element, which sends ultrasonic waves through a large sparse receiver array. The ultrasonic waves are reflect- ed from the imaged target back towards the receiver array, where the signals are recorded and are sent to a processor for further processing and image forming. 29 Figure 2-13: The underwater imaging system with a sparse receiver array (taken from [41]). To achieve the desired imaging resolution of 1mm at the distance of 1m, and with a 3MHz center frequency, the receiving aperture needs to be at least 0.5m x 0.5m. It would of course be impractical to build a full array of this size. The authors decided to build an array with tiles (10 by 10) with smaller sparse sub arrays (0.05m x 0.05m). Each tile contains 32 randomly distributed receiving elements. In the array, there is a total of 3200 elements. The tiled construction of smaller sparse sub arrays is on one hand practical and effi- cient, on the other hand it brings difficulties in precise positioning of the receiving transducers which than greatly affects the image quality. The author developed a method in which the individual receiver elements are cali- brated in the same manner as if they were GPS receivers. A high-precision positioning scanner with a hydrophone attached is used to emulate the role of the GPS satellites in view of the re- ceiver. The hydrophone is used to emit ultrasonic pulse waves, which are intercepted by the receivers. Figure 2-14: Set up for the calibration method with a high/precision hydrophone scanner (taken from [41]). 2.4.2.1 The calibration method If we consider the case, where there are Nh= hydrophone positions, Ne number of ele- ments, then the measured time-of-flight e h from hydrophone position h to element e can be expressed as: τe 𝑕 = te 𝑕 + τe + τ𝑕 + dτe 𝑕 (2.28) where te h is the actual time-of-flight, e and h are the receiver element and hydrophone time- delays respectively (including the transmission and reception channel delays), and dte h is a 30 Gaussian zero-mean measurement noise. It is assumed that e and h are independent of trans- mission and reception angle. Because the same hydrophone is used for all measurements, h is the same for all hydrophone positions. Equation (2.28) is analogical to the pseudorange equation (2.6). The receiver and hy- drophone time-delay components take the place of the clock-error component. The main dif- ference is that (2.28) is expressed in the time-domain, whereas (2.6) is in the spatial domain. The relation between the two is simply the speed of propagation of the waves. Unlike in GPS where the speed of propagation can be assumed constant in all cases (the speed of light), in the underwater imaging application it is an unknown, because the speed of sound differs with temperature (and eventually salinity). That is why the author of this method extended the GPS approach to handle the speed of sound as one of the unknowns. Equation (2.28) can be formulated for each combination of hydrophone position and receiver element number. A system of these equations can be therefore set up and the un- known positions and time-delays can be solved for. Because the actual time-of-flight te h in equation (2.28) is nonlinear, a linearization step is useful, similarly as it is done in the GPS approach. As an initial step, position estimates of the receiving elements (xe0, ye0, ze0) and the speed of sound estimate v0 are chosen, where the actual values are (xe=xe0+xe, ye=ye0+ye, ze=ze0+ze) and (v=v0+dv) respectively. Then the term te h in equation (2.28) can be linearized as follows: te 𝑕 = re 𝑕 v = 1 v (xe − x𝑕)2 + (y e − y𝑕)2 + (ze − z𝑕)2 = 1 v0 + dv (xe0 + Δxe − x𝑕)2 + (y e0 + Δy e − y𝑕)2 + (ze0 + Δze − z𝑕)2 ≈ re0 𝑕 v0 + xe0 − x𝑕 v0re0 𝑕 Δxe + y e0 − y𝑕 v0re0 𝑕 Δy e + ze0 − z𝑕 v0re0 𝑕 Δze + re0 𝑕 v0 2 dv (2.29) where (xe, ye, ze) are the receiver element coordinate error terms, dv is the speed of sound error term, and re0 𝑕 = (xe0 − x𝑕)2 + (y e0 − y𝑕)2 + (ze0 − z𝑕)2 (2.30) is the test distance. The linearized equation set is obtained by expressing the difference between the meas- ured time-of-flight value and the test time-of-flight value: Δte 𝑕 = τe 𝑕 − re0 𝑕 v0 + τe0 𝑕 = xe0 − x𝑕 v0re0 𝑕 Δxe + y e0 − y𝑕 v0re0 𝑕 Δy e + ze0 − z𝑕 v0re0 𝑕 Δze + re0 𝑕 v0 2 dv + Δτe 𝑕 = 1,2,3, . . . N𝑕 , e = 1,2,3, . . . Ne (2.31) where e0 is the time-delay test value for element e, and e=e-e0. By solving the above equation set, a set of position, time-delay, and speed of sound er- ror values are obtained. These can be used to update the previous estimates and the whole process can be iterated until the desired accuracy is reached. This calibration method was perhaps the most inspiring for the author and showed the possibilities of using an established method like GPS being used in a very different environ- ment. Along with the original GPS, it served as a basis for the novel USCT positioning me- thod. The method itself cannot be used for the purposes of the USCT because it makes use of the “high precision scanner” moving the hydrophone to known positions. These known posi- tion values are then used in the equation set (2.31) to calculate the unknown positions of the 31 receiving transducer elements. In USCT however, both the emitting and receiving transducers lie at unknown positions and need to be calibrated. 2.4.3 The RTS calibration approach This calibration approach was developed by Waag and Fedwa for the Ring Transducer System for Medical Ultrasound Research [73] (overview can be found in chapter 2.1). The method takes advantage of the fact that the transducer elements can be used both as an emitter and a receiver. The calibration is also used to determine the sensitivity of individual elements, and so the main concept of data acquisition has already been described in chapter 2.3. A time-shift correction for each element is determined by cross correlation of each rec- orded waveform with a reference waveform 𝑠 (𝑡) (2.4) found by averaging all of the time- aligned waveforms within the range of acceptable sensitivity. Before this step a sinusoidal time-shift bias, caused by the fact that the metal wire is not positioned exactly in the middle of the ring transducer, has to be removed. Figure 2-15: Pulse-echo data before and after calibration. In the plotted data, 13 elements tested low and three elements tested high. The grayscale is linear. (Taken from [73]) Although it is very elegant and lightweight, this method completely omits the geome- trical calibration and only deals with time-delays. The authors assume that the elements on the ring transducer array form a perfect circle. Positioning deviations may actually be partially compensated for by this technique but the calibration results would only be valid for signals coming from the area around the calibration wire. The second and more important reason why not to use this calibration technique for the USCT system in Karlsruhe is the fact that each transducer element must provide the functio- nality to both send and receive. This is not the case in the current Karlsruhe USCT. 2.4.4 Multidimensional scaling Multidimensional scaling (MDS) is a field of study concerned with embedding a set of points in a low-dimensional space so that the distances between the points resemble as closely 32 as possible a given set of dissimilarities between objects that the points represent [14]. It has been popularly used to analyze experimental data in physical, biological, and behavioral sciences. Metric multidimensional scaling (MDS) transforms a distance matrix into a set of coordinates such that the (Euclidean) distances derived from these coordinates approximate as well as possible the original distances. The basic idea of MDS is to transform the distance ma- trix into a cross-product matrix and then to find its eigen-decomposition which gives a prin- cipal component analysis (PCA). Like PCA, MDS can be used with supplementary or illustra- tive elements, which are projected onto the dimensions after they have been computed [24]. As an example, Figure 2-16 shows a map of the United States computed from the driv- ing distances between major cities published in a road atlas, using MDS. Even though the driv- ing distance between two cities is a poor approximation to their actual distance, the resulting map is accurate [7]. Figure 2-16: A map of cities in USA, built by taking the driving distances from a road atlas and applying the MDS algorithm to it. (Taken from [7]) The algorithm is as follows: we first form a square distance matrix D, with the ele- ments di,j being the squares of distance between the i-th and j-th node, di,j = δ2 i,j = (x i − x j) T(x i − x j) (2.32) We then can compute the inner product matrix B: B = − 1 2 JDJ J = I − 1 n E (2.33) where I is the unit matrix, E is the matrix of all ones, and n is the number of the nodes. Now B = XX T , and because B is symmetric positive semi-definite it can be decomposed in to B = VLV T , where L is the diagonal matrix of eigenvalues. We then can compute the matrix of coordinates X = VL 1/2 . The advantage of this method is that it yields a unique solution, where the geometrical axes are set so that the variation of the transducers‟ positions is biggest along the first dimen- sion, the second biggest variation is in the second dimension, and so on. We may then translate and rotate the coordinate system to the desired origin and orientation. This technique could be used to calibrate (find) the positions of transducers in an ultra- sound computed tomography system. The distances between the transducers could be com- puted from the time-of-flight of the ultrasonic signals. The main problem, when applying this technique for position calibration is usually missing data. In order to run the algorithm, the matrix of distances D has to be entirely filled. 33 In the case of the USCT in Forschungszentrum Karlsruhe, the ultrasonic transducer elements are hardwired to operate either as emitters or receivers. They cannot change the operation mode. Thus, the emitter-to-emitter and receiver-to-receiver distances are never recorded. The problem of missing data is a common issue in different technical areas where the MDS algorithm is being used. It has been successfully solved in areas of molecular design using the formalism of distance geometry [1][38][71]. Unfortunately, the molecular design problems are very different in nature from the USCT problem in that only a certain combina- tion of atom constellations are possible to form a molecule. This limitation is exploited in the matrix completion problem. Using a combinatorial approach, all physically impossible con- stellations are avoided, and only the few feasible constellations yield the missing distances among atoms. A different approach to the matrix completion problem was developed for an applica- tion in the microphone array calibration field [7][9]. This method requires building a basis of nodes (microphones), where all inter-node distances can be measured. The number of nodes has to be at least p+1, where p is the number of dimensions of the resulting space (usually p = 3). All other nodes only require having distance measurements available to all of the basis nodes. Other missing distances can be calculated. Unfortunately this approach also cannot be used for the case of the Karlsruhe USCT, because the required basis of at least 4 nodes (transducers) cannot be formed (we are not able to obtain the emitter-emitter and receiver-receiver distances from the time-of-flight of the ul- trasonic signals). 2.4.5 Localization methods in wireless networks Wireless sensor networks got a lot of attention in the past few years. They hold the promise of new applications in the area of monitoring and control like target tracking, intru- sion detection, wildlife habitat monitoring, climate control, and others. Among other problems to be solved in this area of study, the localization of the autonomous nodes (self-organization) is an active area of research [11][39][50]. Although many approaches have been developed, in principle they always try to take advantage of having an infrastructure of so-called beacon or anchor nodes with known loca- tions (and possibly stronger signals) and computing the locations of other nodes relative to these beacon nodes. In some approaches a few nodes are chosen as the anchors in the first phase of the self-localization. Their position is calculated (using Multidimensional scaling or a GPS-like algorithm) relative to other anchors and a base coordinate system is established. The rest of the nodes calculate their own position relative to this coordinate system using only the communication with the anchor nodes. Compared to other applications, the localization methods used in wireless networks have to satisfy requirements on energy efficiency (little computation and communication) and on the other hand do not require such a high accuracy. The main premise of these methods is that the nodes are able to communicate between each other and therefore have both the sending and receiving capability. Moreover the nodes are able to broadcast coded messages with their location. This is not the case in the USCT and therefore these methods are not suitable for the USCT calibration. 34 3 Aims of the dissertation The field of ultrasonic computed tomography has been explored by various groups of scientists for several decades now. Nevertheless, none of these groups has been able to build a fully functional system, which is able to produce repeatable images in non-laboratory condi- tions, which is fast enough for real hospital environment, and which is able to deliver high- contrast and high-resolution images usable in practice. The author had the opportunity to be involved in the Karlsruhe USCT project (de- scribed in detail in chapter 2.2). Also this project presents several unsolved problems. As the title of the dissertation suggests, the primary aim of this dissertation is to develop new calibra- tion methods for an ultrasonic tomography system. The author intends to contribute to the fol- lowing areas: 1. Reconstruction of attenuation images – although the attenuation image reconstruction is quite a broad area (and not specifically linked with the main topic - calibration), the inten- tion is to use this area as a starting point, get practical experience with the project, and ex- plore the possibilities of improving existing attenuation reconstruction methods developed by project colleagues. The individual goals are: a. Get accustomed with the used hardware, software, and methods b. Make measurements on both the 2D and 3D USCT systems in Karlsruhe using cus- tom-built ultrasonic phantoms c. Improve the existing methods of attenuation reconstruction 2. Sensitivity calibration – the currently used image reconstruction methods do not account for sensitivity and directivity differences of the ultrasonic transducers. It is suspected, that the variations in sensitivity may have a big impact on the quality of the reconstructed im- ages. The intention is to create a method for an easy and repeatable calculation of the sen- sitivity parameters of the Karlsruhe USCT transducers, which could then be used (by other team members) in the existing reconstruction algorithms to compensate for the variations. An “empty measurement” (full scan of the s