_: . , :2: : raters This is to certify that the thesis entitled THREE-DIMENSIONAL RECONSTRUCTIONS FROM TWO-DIMENSIONAL BACKSCATTERED ULTRASOUND DATA presented by David Ayres Gift has been accepted towards fulfillment of the requirements for M.S. degree in Computer Science Mdmummumn Date $lq/ (60 0-7639 (fl-‘\\\\ 5 '1“: '7“‘"!iI//”. OVERDUE FINES: 25¢ per day per item RETURNING LIBRARY MATERIALS: Place in book return to remove charge from circulation records THREE-DIMENSIONAL RECONSTRUCTIONS FROM TWO-DIMENSIONAL BACKSCATTERED ULTRASOUND DATA By David Ayres Gift A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Computer Science 1980 ABSTRACT THREE-DIMENSIONAL RECONSTRUCTIONS FROM TWO—DIMENSIONAL BACKSCATTERED ULTRASOUND DATA By David Ayres Gift This thesis describes and analyzes a digital method for the reconstruction of a three-dimensional image from two-dimensional ultrasonic data. The data on which the reconstruction method operates are obtained for the purpose of imaging structural anoma- lies discovered during the nondestructive testing of large objects. It is assumed that the large size of the object being tested con- strains the measurement device to sampling the volume-of-interest from a very limited range of angles. The reconstruction technique is derived from the problem geometry and is based on coordinate system transformations. It is implemented in a modular computer software package which may be easily modified for changes in problem geometry and post-processing requirements, and which allows immediate incorporation of newly- acquired data into the existing reconstruction. The performance of the reconstruction method is tested on a set of data obtained from a phantom of simple structure. ACKNOWLEDGMENTS I would like to thank Dr. Anil K. Jain for the valuable guidance and editorial assistance which he has provided in his capacity as thesis advisor. I would also like to thank the members of the thesis committee, Dr. John J. Forsyth and Dr. Donnie K. Reinhard, for their time, interest and assistance. I wish to acknowledge Dr. John Flora, of Babcock and Wilcox, Inc., for his efforts in formulating the initial problem description. Thanks are also due Dr. E. James Potchen and the Department of Radiology for the use of computing facilities and office equipment, and for the use of data from the Ultrasound Research Project directed by Dr. Reinhard. Finally, I wish to express warm appreciation for the patience, support and assistance of my wife, Debra. ii TABLE OF CONTENTS SECTION Page 1. INTRODUCTION OOOOOOOOOOOOOOOOOO..0...OOOOOOOIOOOOOOOOOOOOOOO 1 1.1 Problem statement ............... ........ ............. 2 1.2 Organization of the thesis ........................... 8 2 O ULTRASONIC BACKS CAHERI NC 0 O O O O O I O C C O O O I C O O O O O O O O O O O C O I I O O O O 1 0 3 O RECONSTRUCTION METHOD 0 C O O O I O O O O O O O O O O O O O O O O O O C O O O O O O O O C O O O O 1 7 O coordinate trans formation 0 C 0 O O C O C O O O C O O C O O O O C O O O O O O O O 19 3 1 3.2 Modification of RS voxel values ...................... 22 3.3 The reconstruction algorithm ......................... 26 4. ANALYSIS OF RECONSTRUCTION METHOD .......................... 32 4.1 Computational complexity and speed of execution ...... 32 4.2 Peint Spread function O...0.00.00.0000...00.00.0000... 38 5. PERFORMANCE OF RECONSTRUCTION METHOD ..... ......... . ........ 45 51 Test data .00...00....OOOOOOOCOOOOOOOOOOOIOOOO0..0.... 45 5 2 Maximumrreplacement rule reconstruction .............. 51 5.3 Summation rule reconstruction ........................ 56 S 4 Sum of squares rule reconstruction ................... 57 5 5 General comments on performance ...................... 59 6. SWARYANDCONCLUSIONS 00......OOOOOOOOOOOOOOOOOOOOOOOOOOO. 62 APPENDIX 0.0.0.000...0..I...OOOOOOOOOOOOOOOOOOOOOOO00.00.0000... 66 REFERENCES OOOOOOOOOOOOOOOIOOOOOOOOOOOOOOOOOO0.00000000000000000 76 iii TABLE 4-2 4-3 4-7 4-8 LIST OF TABLES Page Angles between MS and RS coordinate axes for five views .. 23 Number of operations by class and subprogram of ELECTRA reconstruction program .0.........OOOOOOOOIOOOOO......O 33 Total number of operations for a reconstruction involving N slabs sampled at M points each ............ 33 Resolution—matched PSF for point at (7,8,8); MRR reconstrUCtj-on ......OOOOOOOOCOOOOOOOO......OOOOOOOOOOO 39 Poor-longitudinal—resolution PSF for point at (7,8,8); MRRreconStructionooooooooooooooooooooooooooooo000000o 40 Resolution-matched PSF for point at (7,8,8); SMR reconstruction ..........OOOICOCCOOOOOOO0.000.000.0000. 41 Poor-longitudinal-resolution PSF for point at (7,8,8); SMR reconstruction ......C.............00.0.0.0... ..... 42 Resolution-matched PSF for point at (7,8,8); SSR reconstruction ...0.0.0....00............IOOOOOOOOOOOOO 43 Poor-longitudinal-resolution PSF for point at (7,8,8); SSR reconstruction .....0.........OOOCCOCCOOOCCOCOOO... 44 Vertical view test data .................................. 49 Binarydwindowed vertical view test data; threshold = 1 ... 50 Left view test data ...................................... 52 Binary-windowed left view test data; threshold = 1 ....... 52 Right view test data ..................................... 53 Binarydwindowed right view test data; threshold = 1 ...... 53 Maximum replacement rule test reconstruction .... ......... 54 Binarydwindowed MRR test reconstruction; threshold = 5 ... 54 iv TABLE 5-9 5-10 5-11 5-12 5-13 5-14 5-15 5-16 5-17 Binary-windowed MRR Binary-windowed MRR Summation rule test Binarydwindowed SMR Binary-windowed SMR Sum of squares rule Binary-windowed SSR Binary-windowed SSR Binary=windowed SSR test reconstruction; threshold test reconstruction; reconstruction ..... test test test test test test reconstruction; reconstruction; reconstruction reconstruction; reconstruction; reconstruction; threshold threshold threshold threshold threshold threshold 7 . 10 15 .. 19 75 115 170 Page 55 55 56 60 61 FIGURE LIST OF FIGURES Page Four slabs measured by aiming transducer perpendicular to available surface of test object .................. 4 Three slabs measured from oblique view .................. 6 Scanning views with respect to region of interest ....... 7 Shadowing of reflector R2 by R1 ......................... 14 Measurement raster scanning convention for reconstruc- tion prOblem 0.0.0.0..........OOOOOIOOOOOOOOCO......O. 18 Measurement space (MS) - reconstruction space (RS) coordinate systems for oblique view .................. 20 Reconstruction algorithm flowchart ...................... 27 Voxel scanning order convention (4 x 4 example) ......... 29 Schematic drawing of test reconstruction problem ........ 48 vi 1. INTRODUCTION This thesis addresses the problem of reconstructing a three- dimensional image from two-dimensional data obtained using back— scattered ultrasound. The "image" is actually an approximate map of a function involving the gradients of physical attributes of the interrogated material. The ultrasound echo data are obtained in a constrained measurement geometry: the region of interest lies within a piece of material of which only one large flat, or slightly curved, surface is accessible to the measurement apparatus. Thus, the echo data represent backscattered ultrasound which "sees" the region of interest from only a small range of angles. The problem just described is typical of applications in non- destructive testing (NDT) of materials [1,2]. This testing is some— times done using X-rays, but is more often carried out with ultra- sound if the material being tested is a metal or plastic which is highly absorptive for X-rays. Also, it is often the case that the test object is a part of a large structure or assembly which allows only limited access to the measurement apparatus (in this case, an ultrasound transducer), thus precluding the use of transmission techniques. The goal of ultrasonic nondestructive testing is to detect and image cracks, stress defects, and other flaws in the material, as well as detecting and imaging "intentional" structures such as holes and fasteners. Ultrasound is also gaining increasing popularity as a medical imaging tool because it avoids the use of ionizing radiation and has the potential to provide information on the histological structure of the imaged tissue. Much of the recent theoretical work on ultra- sonic scattering has developed from medical applications [3,4]. While the primary subject of this thesis is materials testing, occasional reference will be made to the clinical use of ultrasound for the purpose of comparison. The reconstruction technique proposed herein is based on a spatial coordinate system transformation. The method requires no special assumptions regarding the physical nature of ultrasonic reflectivity and may be implemented in highly modular computer soft- ware, making it adaptable to changes in measurement geometry and post-processing requirements. 1.1 Problem statement. Within a region of interest (ROI) contained in a field on mate— rial, one desires to image artifacts which deviate from the homo- genious nature of the surrounding material. The desired image, furthermore, is to be three-dimensional -- allowing the observer the capability of viewing the inhomogeneities in the ROI from any direc- tion, such that their individual shapes and extents, and their posi- tional relationships to each other may be ascertained. Interrogation of the ROI is accomplished using ultrasound. In this method, the ROI is "illuminated" by a narrow-bandwidth ultra- sound beam. Wavefronts reflected from within the ROI are measured with a transducer, raster-scanned across a "scanning aperture". The signal at the aperture is sampled discretely by the raster- scanned transducer. The image obtained from the echo data is, therefore, specified by a two—dimensional matrix of discrete picture elements (pixels). In the case of the measurement apparatus under consideration here, however, the output of the scanned transducer is time-gated and integrated over the duration of the gate. Thus, each sampled point in the echo image results from a train of reflections returning from a region of finite thickness, rather than from a plane. The thickness of the region actually measured is equal to the product of the duration of the time gate and the speed of ultrasound in the interrogated material. This thickness represents the longitudinal spatial resolution of the measurement apparatus. The lateral extent of the region measured is determined by the extent of raster scanning. The lateral resolution of the measurement apparatus is equal to the lateral resolution achieved by the interro- gating ultrasound beam or by the raster unit of the sampling device, whichever is the larger. This measurement environment is shown schematically in Figure 1-1. The measured data derived from the sampled images take the form of "slabs" passing through the ROI. The face of each slab is divided into a matrix of pixels - the discretely sampled points of the echo image; each pixel has a finite depth -— the longitudinal resolution, or slab thickness. Because of the spatial depth of the image matrix elements, they are more properly thought of as volume elements, or "voxels", and shall be referred to as such henceforth. Transducer scanned here 7 L A I L A ' fir / Figure 1-1. / Four slabs measured by aiming transducer perpendicular to available surface of test object. Measurement of the ROI is subject to the constraint that the ROI is contained within a field of material, only one face of which is accessible to the measurement apparatus. Thus, the ROI may be viewed from only a limited range of angles. One such oblique view is shown in Figure 1-2; the top face of the cube containing the ROI is the surface accessible to the measurement device. As depicted in Figure 1-2, a wedge coupled to the accessible surgace with gel may be used to establish the appropriate angle for an oblique view. The choice of the angle at which an oblique view is made is arbitrary, but subject to some constraints. The choice is limited by the extent of the available surface of the test object. A larger surface area allows views at greater (more oblique) angles than does a smaller surface. Ideally, one would like to observe the ROI from all sides, and large oblique angles are better in this sense than are small oblique angles. However, the area where the wedge is coupled to the available surface is a reflecting interface. As the angle of the view increases, so does the amount of reflection and refraction of the interrogating and echo signals at that interface. In addition, more oblique views place the ROI at a greater distance from the trans- ducer than do less oblique views and the ability of ultrasound to penetrate the test material becomes a critical factor. While the ROI may be viewed from any point on the available sur- face, we shall consider the hypothetical case in which five views of the ROI are obtained as in Figure 1-3. Here, the ROI is measured from a point on the surface where the line from the center of the R01 to the center of the scanning aperture is orthogonal to the plane of the surface -- a "vertical" view, designated "V" in this figure -- and Transducer 1. scanned here 1/ / / Figure 1-2. / l l l l )— Three slabs measured from oblique view. interest Figure 1-3. Scanning views with respect to region of interest. from four points offset from the vertical view along the arms of a cross centered at the vertical view position, such that the ROI is observed on four "sides". We shall refer to these views as "front" (F), "left" (L), "back" (B), and "right" (R). The "front" direction is arbitrarily chosen. The reconstruction problem may now be simply stated: Measure- ments are made of the properties of a region of interest within a field of material. The measurement data represent two-dimensional slabs passing through the ROI. The slabs are each comprised of a matrix of pixel values (lying on the face of the slab); each pixel value is the value of the property of the ROI as measured at the pixel location. The longitudinal resolution of the measurement appa- ratus is finite, so the pixels, and the slabs, have finite depth or thickness. Measured slabs are obtained at several angles through the ROI. From these data, a three-dimensional reconstruction of the ROI is desired. 1.2 Organization of the thesis. This thesis is organized as follows: Section 2 briefly discusses the nature of ultrasonic backscattering in the materials testing con- text. Section 3 describes the pr0posed reconstruction method and its implementation as a computer program. The reconstruction method is analyzed in Section 4, with attention being given to its computational complexity, speed of execution, and an empirical assessment of the point spread function. The results of testing the performance of the reconstruction method on data obtained from a phantom of simple 9 structure are described in Section 5. Section 6 summarizes the contents of the thesis, discusses the conclusions to be drawn from this work, and suggests tOpics for further investigation and development of this reconstruction method. 2. ULTRASONIC BACKSCATTERING Before proceeding with a detailed description of the geometry of this reconstruction problem and the proposed solution, it is important to consider the physical prOperties which are measured by backscat— tered ultrasound, as well as the nature of ultrasonic scattering. Ultrasound consists of a train of pressure wavefronts propagating radially outward from their source. These waves may be focussed, collimated, refracted, diffracted, and reflected by the properties of the media through which they travel. We shall be concerned primarily with backscattering of ultrasound. As interrogating ultrasound waves travel away from their source and encounter a reflector, they are scattered by the reflector in many directions, depending on the orientation and other characteristics of the reflector. Those waves which are scattered back toward the source, in a direction rotated 180 degrees to the direction of propagation of the interrogating waves, constitute the "backscattered" ultrasound. In this section we consider the following questions: What con- stitutes a reflector, ie. what physical properties of the transmit- ting medium cause scattering? What is an "ideal scatterer" and what is the nature of scattering from it? What deviations might be expec- ted in the case of non-ideal scatterers? Gore and Leeman [3] demonstrate that, if one measures backscat- tered waves only, the nature of the echo is dependent on the shape and 10 11 spectral content of the interrogating pulse, and also on the differ- ence of compressibility and density gradients in the interrogated medium. In addition, since the incident pulse is of finite longitudi- nal and transverse extent, the measured echo at any point in time is actually due to the average of backscattering over the volume defined by the extent of the incident pulse at the time of the genesis of the echoes. If the echo is time gated as it is being sensed, there is further smoothing of the recorded echo image due to finite time reso— lution. Thus, although ultrasonic backscatter is caused by density and compressibility gradients within a medium, one may only measure the pulse and time smoothed difference of such gradients. Neglecting any time constraints due to sampling of the echo, two scatterers hav- ing different compressibility and density gradients but the same dif— ference between these gradients may produce equivalent echo amplitudes. Recently, Aks [3] has suggested that viscosity gradients also produce ultrasonic scattering. We may ignore this source of scatter- ing in the problem at hand, since we are considering defects in the structure of a metal field. Two properties of ultrasonic backscattering are important with respect to the imaging problem at hand: (1) anisotropy of backscat- tering from a non-ideal scatterer; (2) non-additivity of backscatter- ing from multiple longitudinally-colinear scatterers. In the first case, consider a reflector embedded in a non-disper- sive homogeneous medium. A transducer is used to beam interrogating ultrasound toward the reflector, and to sample the backscattered sig- nal. The transducer is also embedded in the ideally—propagating medium, but is restricted to be located only at the locus of points 12 equidistant to the reflector. If the backscattered signal produced by the interrogation of the reflector is the same for any position of the transducer, backscattering from the reflector is isotropic. A reflector for which backscattering of ultrasound is isotropic is labelled a "Rayleigh scatterer" and is said to produce monopole scattering. A scatter will appear as a Rayleigh scatterer if and only if two conditions are satisfied: (1) there are no viscosity, com- pressibility, nor density variations within the scatterer; (2) the size of the scatterer is small with respect to the wavelength of the inter- rogating ultrasound [5]. Even under the assumption that the first condition is always satisfied, isotropy of backscattering may not gen- erally be expected to occur in a materials testing context. As an example, the frequency of the interrogating ultrasound is typically 2 MHz [6]. The speed of ultrasound in steel is 6900 meters/second [5]. Thus, the wavelength typically encountered in nondestructive testing for structures in steel is 2.95 mm. For a scatterer to appear as a Rayleigh scatterer its largest extent must be much smaller than 2.95 mm. Structures violating this size limit are quite easy to imagine. If a structure —- a crack, for example —- is several times larger than this limit, its scattering becomes very "directional". A crack which is long (with respect to the wavelength of the ultrasound) and which is oriented orthogonally to the direction of propagation of the interrogating signal will produce a larger backscattered signal than will an identically "long" crack oriented obliquely to the direction of propagation. In the extreme case, a long crack seen "end-on" by the transducer will backscatter relatively little of the interrogating 13 beam. The key problem raised by the inability to assume backscattering isotropy for all reflectors in the ROI is this: What does the absence of a reflected signal mean? Does it indicate the absence of reflec- tors? Not always -- it may simply fail to indicate the presence of reflectors which are oriented obliquely to the transducer's direction of aim. Indeed, a large obliquely-oriented crack could return the same (or possible a lesser) echo signal as would a small orthogonally- oriented crack. Thus, a large echo signal indicates a strong reflec- tor, but the information carried by a weak or zero echo signal is ambiguous. The second property of ultrasonic backscattering which is impor- tant to consider here involves the nature of the echo signal received from multiple reflectors which are colinear with each other and with the transducer. For the moment, consider two such reflectors (R1 and R2) embedded in ideally-propagating (ie. lossless) material as shown in Figure 2-1. A transducer is also embedded in the material and is aimed so that R1 is directly in front of R2. The transducer fires an interrogating pulse, then receives the echoes returned from R1 and R2. The echo from R1 will be strong, but the echo from R2 will be weak, even if R1 and R2 have equal reflection properties. This is due to two things. First, the strength of the interrogating pulse reaching R2 is reduced from that reaching R1 by the amount of the signal scat- tered by R1. In other words, R1 "shadows" R2. Second, the echo returned from R2 will also be attenuated when it traverses R1 as it travels toward the transducer. R1 will scatter a portion of this echo signal, also. 14 Transducer * 1 / R2 R1 Figure 2-1. Shadowing of reflector R2 by R1. 15 We can call this second characteristic of ultrasonic scattering "shadowing". Two identical reflectors, oriented identically with res- pect to the transducer will not produce identical reflections if one reflector is shadowed by another. Indeed, if the scattering surface of the shadowing reflector is obliquely oriented to the direction of propagation of the interrogating beam, little or no interrogating ultrasound may reach the second reflector at all (if the oblique, shadowing reflector is large with respect to the interrogating wave- length). Thus, anisotropy and shadowing can combine to hide reflec- tors by causing attenuation of the interrogating beam, the reflected signal, or both. This discussion has ignored the interference phenomena which are due to the wave nature of ultrasound. It has been mentioned that the interrogating ultrasound in this problem is pulsed. The pulses are of finite duration and will, therefore, result in pulsed echoes of finite duration. If echo pulses from two reflectors overlap in time as they reach the transducer, their wavefronts will constructively or destruc- tively interfere to produce either a larger or smaller amplitude sig- nal. Two reflectors lying nearly equidistant from the transducer and within the transverse extent of the interrogating pulse could be seen by the transducer as either a single strong reflector or a single weak reflector, depending on the size of their separation with respect to the interrogating wavelength. To summarize, in an environment where non-ideal reflection may be expected to occur, weak echoes (or the absence of echoes) are ambigu- ous as indicators of the presence and characteristics of reflectors. Strong echoes, on the other hand, indicate the presence of a reflector 16 or, at least, the presence of small multiple neighboring reflectors whose echo pulses are (1) not resolved by the measurement apparatus, and/or (2) overlap to constructively interfere. 3. RECONSTRUCTION METHOD Reconstruction of three-dimensional images from ultrasound data has been pursued in several ways. Holography is perhaps the most elegant approach, utilizing phase information ignored by other methods [7-10]. The reconstruction of holographs by numerical means is time and memory intensive, however, and optical reconstruction requires complicated measurement apparatus design. Amplitude-based methods can make more direct use of the measured data. Norton and Linzer [11,12] have preposed two systems of transducers -- one circu- lar and one spherical -- for the production of two and three-dimen- sional images. These designs rely on the ability to fit the region of interest inside the measurement device. The method of solution proposed for this reconstruction problem under constrained measurement geometry is derived from the problem geometry, and is essentially equivalent to a technique described by Onoe, et al [13], for the display of reconstructed images. Figure 3-1 presents, once again, the sampling scheme used to obtain echo data. A transducer is raster-scanned across the face of the object being investigated. Data are time-gated such that a completely-scanned rectangle on the surface of the object yields reflection data for a volume slab which intersects the region of interest. In this case, we assume that the transducer is scanned across an offset such that the slab's major face lies orthogonal to the direction of propagation 17 18 Scan lines #1 #2 8 Measured slab \ \ Figure 3-1. Measurement raster scanning convention for reconstruction problem. 19 of the interrogating ultrasound (ie. not necessarily parallel to the available surface of the material). Echo data are sampled at discrete intervals along each scan line, so the scanning geometry is well- suited to description by cartesian coordinates. We shall refer to the coordinate system describing each slab of sampled data as measurement space (MS), and denote the MS coordinates as xm, ym, and zm. Figure 3-2 illustrates the relationship of measurement space to the region of interest. We may arbitrarily describe the region of interest using cartesian coordinates, and will henceforth refer to this coordinate system as the reconstruction space (RS), denoting the RS coordinates as x, y, and 2. From Figure 3-2 we see that any set of MS coordinates may be transformed into RS coordinates by translation and rotation. Further— more, the translation and rotation may be decomposed into translations along, and rotations about, each of the RS axes. The reconstruction problem may now be treated in two parts: (1) transformation of the MS coordinate values for each measured voxel into the RS coordinate values; (2) modification of the RS voxel data value by the newly- measured data. The next sections discuss each of these reconstruction steps. 3.1 Coordinate transformation. Translation-rotation transformations may be carried out by a matrix multiplication involving direction cosines [14]. Referring to Figure 3-2, let the origin of the MS coordinate system lie at RS coordinates (xo,yo,zo). Notice that 20 §_0. The rotation of MS axes 20 (XO,YQ,ZO) Figure 3-2. Measurement space (MS) - reconstruction space (RS) coordinate systems for an oblique view. 21 with respect to RS axes is given by the direction of offset of the scanned aperture from the region of interest (V, F, L, B, R -- see Figure 1-3) and the angle A between the direction of transducer aim and the z axis of RS. Specification of the rotation may be decomposed to the angles between each MS-RS coordinate axis pair. Letting Aab be the cosine of the angle between axis a and axis b, a matrix of these "direction cosines" may be denoted as A A xxm xym xzm Letting M = (xm,ym,zm)T 3 (X9Y9Z)T ; 2'1 II T T = (xo.yo.zo) ; the complete coordinate transformation may be written R = DM + T . or, augmenting D and M such that C II a [DIT] a T and Ma (Xm,Ym.Zm.1) 9 then, R D M . 22 Table 3-1 shows the inter-axis angles for the five views illus- trated in Figure 1-3. As an example of the direction cosine matrix used for the coordinate transformation, consider a "right" (R) view for which the transducer is aimed at 30 degrees with respect to the z axis of RS, and the center of the scanning aperture lies 10 units along the y axis and -17 units along the z axis (see Figure 3-2): a _cos(90) cos (30) cos (90+30)_ D = cos(180) cos(90) cos(90) cos(90) cos(90-30) cos(30) ‘- = -1 o o o .5 .366_ T = (0,10,-17)T and To .866 -.5 0" D =-1 o o 10 a o .5 .866 -17__ 3.2 Modification of RS voxel values. Having transformed the coordinate values from MS to RS, one must now modify the existing data for the appropriate RS voxel in light of the newlydmeasured data obtained for the transformed MS voxel. As discussed in Section 2.2, a large echo signal is a more reliable 23 Table 3-1. Angles between MS and RS coordinate axes for five views. Angles between MS & RS axes (degrees)_ VIEW ANGLE Km ym zm RS V 00 0 90 90 x 90 0 90 y 90 90 0 2 F A; 0 9o 90 x 90 A 9O+AF y 9 0 9 0- F AF Z B A3 180 90 90 x 90 180-AB 9O—AB y 90 9o-AB AB 2 O 180 90 90 y 90 90-AR AR 2 L A0 90 180-A 90- x L L o 90 9oAL y 90 9o-AL AL 2 24 indicator of the presence of a reflector than is a small echo signal. On the basis of this statement, one might prOpose a very simple voxel- value modification rule: New voxel value = maximum(current value, newly-measured value) . We shall call this rule the maximum replacement rule, or MRR. Due to suboptimal spatial resolution in the measurement apparatus and to computational errors (round-off, etc.) within the reconstruc- tion algorithm, the image of a point reflector located at a specific voxel in RS may spread into neighboring voxels. The Spreading of the image of a point reflector is described by the point spread function (PSF) of the reconstruction technique. A potential flaw of the MRR is the granting of equal "reflector-presence credibility" to voxels affected by the PSF of a true reflector as it does to voxels in which a true reflector lies. The MRR ascribes an equal RS voxel value to a voxel for which the values 0, O, O, O, 10 have been measured (from five views) as it does for a voxel for which the values 10, 10, 10, 10, 10 have been measured. Thus, some scheme for taking into account all observed values for a given voxel seems desirable. A summation rule (SMR) in which the RS voxel value V is equal to the sum (over all measurements) of the values observed for that loca- tion, vi, V = Xvi , has the advantage of reinforcing values for multiple observations of a reflection at the same location, and requires only as much memory as is needed for the reconstruction array. 25 The SMR weights small and zero measured values equally with large measured values, and this property runs counter to the nature of ultrasonic scattering. Large observed values should be allowed to have greater influence on the reconstruction than small or zero values due to the less ambiguous interpretation of reflector presence offered by a large signal than a small signal. Weighting larger signals more heavily than smaller signals may be accomplished by a sum of squares modification rule (SSR), 2 V- 2V1 o The SSR weights each observed value by itself and, like the SMR, requires only as much memory as is needed to store the reconstruction array. Voxel-value modification rules which take the number of observa- tions, n, of each RS voxel into account and are similar to the SMR and SSR are the arithmetic mean V = (1/n) Xvi and mean square V = (1/n) Xvi rules, respectively. These rules require twice as much memory as is needed for the reconstruction array, since the number of observations of each R8 voxel must also be stored. A reasonable reconstruction volume may be 128 x 128 x 128 voxels, or 2000k values (2M values). Computers with 32 bit word size are commercially available at moderate cost. These machines are, in general, designed to address 2-4M words 26 of main memory. Mean-based rules would require 4M words of storage for data alone for a 1283 RS array. The cost of the extra memory requirement for mean-based rules is prohibitive - either in dollars for main memory, or in time for the use of disc memory and the additional software overhead of a scheme to swap the RS data in and out of main memory. For this reason, we shall restrict our consideration of voxel value modification rules to those requiring a minimum of memory, namely, maximum replacement (MRR), summation (SMR), and sum of squares (SSR). 3.3 The reconstruction algorithm. As previously discussed, each slab of measured data may be des- cribed with respect to a coordinate system whose z axis passes through the slab's center and is orthogonal to the major plane of the slab. In this case, the slab may be transformed, voxel by voxel, from MS to RS, where the chosen voxel-value modification rule is applied. A flowchart of the reconstruction algorithm appears in Figure 3-3. Each step of this algorithm is easily implemented as a subroutine, allowing the algorithm as a whole to conform to changes in data col- lection, transformation, reconstruction-voxeldvalue modification rule, and post-processing of the reconstructed image. Following the flowchart in Figure 3-3, the algorithm begins by initializing the elements of the reconstruction array to a value selected to denote "no data" (eg. -1). Next, the data for one measur- ed slab are read in. 27 c ) Initialize reconstruction array elements to -1 3L 'I . Read next slab data J Compute (or look up) augmented matrix of direction cosines for current slab. Generate next voxel coordinatej (xm.ym.zm) Transform coordinates from MS I to RS (x,y,z) iApply selected modification J rule to RS voxel value More voxels in yes current slab? no ‘ other slab yes available? no Post-processing & output I of reconstruction array END Figure 3-3. Reconstruction algorithm flowchart. — .—._.__._.__ 28 The input data consist of: a) slab identification number (unnecessary to the reconstruction -- for bookkeeping purposes only); b) the depth of the slab from the transducer (zm -- known from the product of the speed of sound in the interrogated material and the sampling delay divided by two); c) the coordinates of the transducer position in RS (x0, yo, zo —- x0 and y0 are measured on the surface available to the transducer, -zo is the depth of the top of the region of interest from the available surface -- see Figure 3-2); d) the view, or measurement position, code (V, F, B, L, R -- determines the nature of axial rotation specified by the direction cosines used in the coordinate transformation -- see Figure 3-2); e) the angle of transducer aim with respect to the RS 2 axis (the z axis is chosen to be orthogonal to the available surface at the vertical (V) view position); f) measured data values (in the order scanned —- raster scanning is assumed to proceed left-to-right (-xm to +xm) across the scanning lines, the lines are scanned in order from the edge of the scanning aperture nearest the V-view aperture to the edge furthest from the V-view aperture (-ym to +ym) —- see Figures 1-3 and 3-4). Next, one computes or looks up the elements of the matrix of direction cosines appropriate for the transformation of coordinates for the view from which the current slab‘s data were obtained. 29 0 Scanning direction:— COLUMNS: O 1 2 3 ROWS 1 2 3 4 O 5 6 7 8 (l.5,-O.5) 1 O 9 10 11 (.5,0) 12 (1.5,0) 13 14 15 16 ym (zm into page) Figure 3-4. Voxel scanning order convention (4 x 4 example). 30 Following this, a loop is entered and traversed separately for every input measured data value. Within the loop, the number of the current measured voxel value (equal to the index of the loop) is used to compute the voxel's MS coordinates using the assumed scanning order convention (-xm to +xm and -ym to +ym; see Figure 3-4). These MS coordinates then are multiplied by the augmented direction cosine matrix to transform the MS coordinates to RS coordinates. The RS coordinates identify the apprOpriate reconstruction array element whose value is then modified by the measured voxel value in accordance with the selected modification rule. This 100p is repeated until all data values for the current slab are processed into the reconstruc- tion. If data for another slab are available, these data are read and the reconstruction process is repeated. If no additional input data exist, the reconstruction is complete and the reconstruction array is ready for post-processing and display. ELECTRA ("[picture-volume] ELEment Coordinate-Transform Reconstruction Algorithm") is a FORTRAN IV computer program which implements the reconstruction algorithm illustrated in Figure 3-3. A listing of the ELECTRA code may be found in the Appendix. The driver (main program segment) is responsible for managing the flow of the reconstruction algorithm, making branching decisions, and defining basic variables and storage requirements. (For simpli- city, the program is written for reconstruction of a 16 x 16 x 16 matrix, or a 4k—voxel volume.) The driver calls subroutine DIRCOS to compute the matrix of direction cosines apprOpriate for the slab's measurement position and 31 the angle of transducer aim. For each slab voxel in order, subroutine VOXGEN is called with the current voxel number and depth of the slab from the transducer. VOXGEN returns the MS coordinates of the current voxel. Subroutine TRNSFM is then called with these coordinates and the augmented direction cosines matrix. TRNSFM performs the spatial transformation from MS to RS and returns the voxel's RS coordinates. The RS coordinates and the data value for the current voxel are passed by the driver to the user-selected voxel value modification subroutine, which alters the reconstruction array. ELECTRA can provide output of the completed reconstruction array either as tables of absolute values of array elements or as "binary- windowed" tables consisting of two symbols: "@", representing values greater than or equal to a user-specified threshold, and ".", repre- senting values less than the threshold. 32 4. ANALYSIS OF RECONSTRUCTION METHOD 4.1 Computational complexity and speed of execution. The analysis of the computational complexity of ELECTRA is carri- ed out by tabulating the number of computational operations performed in each subprogram for a reconstruction utilizing the maximum- replacement voxel value modification rule (MRR). This tabulation is presented in Table 4-1. The entire table is traversed once for each measured slab. The inner loop of processing routines - VOXGEN, TRNSFM, and MRR -— is traversed once for every slab volume element, or voxel (sampled point in the raster scan). Table 4-2 displays the totals for each type of operation, in the event that N slabs of measured data, each sampled at M points (M = product of the number of scanning lines and the number of raster units per line), are presented for the reconstruction. The N slabs represent the total number of slabs used; each of these slabs may have been obtained at different angles with respect to the region of inter- est. The numbers and expressions shown in Tables 4-1 and 4-2 are exclusive of operations performed for the initialization and output of the reconstruction array. They also represent a worst-case analysis in that certain of the counted operations will not always occur. This latter condition pertains, for example, to some of the operations 33 Table 4-1. Number of operations by class and subprogram of ELECTRA reconstruction program. OPERATION CLASS SR ‘AS' yp_ OTHER For each slab (N slabs): Driver .......................... - - - 1 disc access. DIRCOS .. ..... ........... ...... .. 12 3 1 5 comparisons, 4 COS calls. Driver .......................... 3 - - -- For each voxel (M voxels/slab): VOXGEN ..................... 1 6 3 -- TRNSFM .... ...... ........... 3 9 9 -- Driver .......... ........ ... 3 — - -- MRR ... .......... ........... - 6 — 10 comparisons or logical, 1 AMAXl call. SR = simple replacement, eg. x = y AS = add/subtract MD = multiply/divide Table 4-2. Total number of operations for a reconstruction involving N slabs sampled at M points each. OPERATION CLASS NUMBER OF OPERATIONS SR N((12 + 3) +-M(1 + 3 + 3)) = N(15 + 7M) AS N(3 + M(6 + 9 + 6)) = N(3 + 21M) MD N(l + M(3 + 9)) = N(l + 12M) Logical/ comparison N(5 + 10M) Function calls N(4(COS) + M(AMAX1)) I/O N 34 counted for subroutine MRR, which are not performed if the transformed measurement space (MS) voxel coordinates are determined to lay outside the reconstruction space (RS) boundaries. The operation counts tabulated in Tables 4-1 and 4-2 are grouped into classes of operational types which are assumed to execute in roughly the same amount of time on a typical minicomputer. A simple replacement (SR) operation which copies a value from one memory loca- tion to another via a Load-Register, Store-Register instruction sequence is expected to execute rapidly. Add/subtract (AS) operations are assumed to be equivalent in execution speed (Load—Register, Add or Subtract, Store-Register) and slower than SR Operations, as are multi- ply/divide (MD) Operations (Load-Register, Multiply or Divide, Store- Register). Logical and arithmetic IF operations (Load-Register, Subtract, Jump-If) are assumed to require about the same amount of time as MD operations (and to be slower than AS operations). Operations proceeding at relatively slower rates than any of the above are calls to functions and accesses to external memory devices. Notice that for all classes of Operations the order of computa- tional complexity is NM; M = jk, where j is the number of raster lines per slab scan and k is the number of raster units (sampled points) per raster line. The development and testing of the ELECTRA software has been done on a Varian V75 minicomputer (currently known as the Univac V77-600). This machine has a word size of 16 bits, an effective CPU cycle time of 195 nsec and is equipped with a floating-point processor. To ob- tain an estimate of the time required to perform a given reconstruc- tion, we shall use this machine as a reference processor, under the 35 following assumptions and stated system properties [15]: (1) Assumption: All operations involve single-precision (32 bit) floating-point arithmetic. (2) Property: Average time required (ATR) for SR operations is 4 microseconds (usec). (3) Property: ATR for AS operations is 7 nsec. (4) PrOperty: ATRfor MD operations is 9 nsec. (5) Assumption: ATR for an arithmetic IF or logical comparison is equivalent to the ATR for a MD operation, 9 usec. (6) Assumption: ATR for calls to the functions COS and AMAXl is long, 100 usec. (7) Property: ATR for disc access is 55 msec. Not all Operations performed in ELECTRA are on floating-point quantities. Since integer Operations proceed at a faster rate than floating-point operations, assumption (1) will yield a conservative estimate of the time required for the reconstruction. This will be Offset, however, by the exclusion from consideration of jumps to and from subroutines and functions, time required for the memory manage- ment of large matrices, and operating system overhead. Using these figures and estimates for the ATR for execution of each operation class, we may estimate the time required for a given reconstruction problem: Let each slab be raster scanned, sampled at 128 points on each of 128 scan lines. The number of voxels measured per slab, then is M = 128 x 128 = 16 384 voxels/slab. Suppose, also, that the region of interest is interrogated from 9 36 views (vertical, plus 2 views at different angles each from positions front, back, left, and right of the region's center) and that 10 slabs (or 10 time gates of the backscattered signal) are sampled from each view. Thus, the total number of slabs involved in the reconstruction is N = 9 x 10 = 90 slabs. With these values for N and M, the number of Operations Of each class required to perform the reconstruction is: SR 10 323 270; AS. 30 966 030; MD: 17 694 810; Logical/comparison: 14 746 050; Function calls: 1 474 920; I/O (disc access): 90. Using the ATR values cited above for the reference computer, the estimated average time for a reconstruction from 90 slabs sampled at 128 x 128 points each is approximately 702 seconds, or 11.7 minutes. Test runs of ELECTRA on the V75 for reconstructions involving 104 slabs of 16 x 16 (256) voxels have shown that the actual time required for reconstruction is approximately twice that predicted by the above complexity analysis. The V75 runs under the control of a real-time multitasking executive. Although no other user jobs were running during the test reconstructions, the Operating system suffers a real time clock interrupt every 5 msec, in response to which it services interrupt handlers and peripheral drivers, updates system tables, and re-evaluates next-jOb-to—execute status. The overhead imposed on the reconstruction time by these executive duties cannot be directly 37 measured, but accounts for most of the difference between the observed and predicted performance. In a non-destructive testing environment where reconstruction speed is critical, a dedicated processor would more nearly match the calculated execution speeds. Furthermore, the reconstruction method implemented by ELECTRA proceeds slab-by—slab. If each slab's data were input to ELECTRA as they became available at the output of the measurement apparatus, the reconstruction would be complete shortly after the final slab's data were measured. Processing of the ' "current" slab data during measurement of the "next" slab data would greatly reduce the time required for a reconstruction. Computational complexity of the SMR and SSR rules may be discus- sed with respect to that of the MRR. These rules differ only in the modification rule computation. Thus, where the MRR has one AMAXl function call with two arguments, the SMR has one addition Operation, and the SSR has one multiplication and one addition operation each. Speed of execution for the MRR was computed under the assumptions that (on a reference processor) the average time required (ATR) for AMAXl is 100 msec, the ATR for an addition is 7 usec, and the ATR for a multiply is 9 usec. Thus, the SMR and SSR modification rules are~ substantially faster than the MRR. Using the analysis for a 1283 reconstruction from 90 slabs of 1282 voxel measurements, we find that the SMR.will require 9.4 minutes to achieve the reconstruction, and the SSR 9.6 minutes. These figures are to be compared to 11.7 minutes required by the MRR. 38 4.2 Point spread function. The spreading or blurring of the image of a point object is des- cribed by the point spread function (PSF) of the imaging or image- reconstruction technique. Two conditions account for blurring of a point image as reconstructed by ELECTRA: (1) If the reconstruction array has finer spatial resolution than does the measurement apparatus, one can expect a spreading of the image of a spatially small reflector over several RS voxels. (2) Quantization error produced by rounding— off computed RS coordinates to integer values used in addressing ele- ments of the reconstruction array could, on occasion, cause errors of i1 value along any coordinate direction in RS, resulting in the map- ping of an MS voxel to an adjacent neighbor of the correct RS voxel. The PSF resulting from these sources of error was studied empiri- cally. The slab data for the reconstructions were generated by hand using a graphical method inherently different from the coordinate- transform reconstruction method; all observations of the point reflec- tor were assumed to occur in the far-field of the transducer, where the beam energy is spatially homogeneous. The data set consisted of 104 slabs, each of 16 x 16 voxels, measured from 5 views of the point object -- V, at (xo,yo,zo) = (7.5,7.5,-0.5); F, at (7.5,15.5,-0.5); B, at (7.5,-0.5,-0.5); L, at (-0.5,7.5,-0.5); R, at (15.5,7.5,-0.5). We first consider the PSF due solely to computational error, where the spatial resolution of the MS and RS are equal. We shall call this the "resolution-matched" case. Table 4-3 displays the re- construction Of an ideal point reflector of strength "10" at (x,y,z) = (7,8,8). Under the maximum replacement rule, the point has spread 39 over 3 voxels, and has acquired an "L" shape. Loss of measurement spatial resolution in NDT applications is generally in the longitudinal direction. This is due to the high speed at which sound travels in metals (such as steel) typically en- countered in NDT problems. We consider next the case where lateral resolution in MS and RS is matched, but MS longitudinal resolution is 3 times less than RS longitudinal resolution. That is, a point physi- cally 1 voxel in all dimensions is measured as being 1 voxel-unit wide in the xm and ym directions and 3 voxel-units deep in the zIn direction (the measured slabs are three times "thicker" than the length of a RS voxel side). For the point at RS coordinates (7,8,8), the PSF under these conditions is shown in Table 4-4. Table 4-4 illustrates the PSF caused by the combined effects of measurement error and computational error. Table 4-3 shows the PSF Table 4-3. Resolution-matched PSF for point at (7,8,8); MRR recon- struction. Z-SECTION 8: X: 0 A 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Y/ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 10 10 0 0 0 0 0 0 0 0 9 0 0 0 0 0 0 0 10 0 0 0 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0 0 O 0 0 0 0 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 40 (Y = 12 through 15 all zeroes) Poor-longitudinal-resolution PSF for point at (7,8,8); MRR reconstruction. Table 4-4. Z-SECTION 8: 9 10 11 12 13 14 15 8 2 1 0 x. 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 11 oooooooomooo 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 0112345678901 11 9 10 11 12 13 14 15 8 1 2 0 Z-SECTION 9: 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 oooooooomooo 000000000000 oooooooomooo 000000000000 000000000000 000000000000 000000000000 000000000000 012314567890]. 11 9 10 11 12 13 14 15 8 2 Z-SECTION 10: 0 1 X: 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 oooooooomooo 000000000000 000000000000 1 .I. 1 000000000000 oooooooomooo 000000000000 000000000000 000000000000 000000000000 000000000000 012345678901 11 41 due to computational error alone. It is seen that in either case, an ideal point reflector is imaged as a complicated configuration of points resembling a lopsided asterisk. The maximum replacement rule gives equal image credibility to all of the voxels to which the image of the actual point has spread. The effect of this PSF on images of irregularly-shaped reflecting objects could render those images mean— ingless. The PSF of the summation rule under resolution-matched conditions appears in Table 4-5. The PSF under conditions of poor longitudinal resolution appears in Table 4-6. The actual location of the point reflector in each case is at (7,8,8); the point's echo has a measured value of 10 from each View. The spatial shape of the PSF for the summation rule is identical to that of the PSF for the MRR. The PSF of the summation rule is not Table 4-5. Resolution-matched PSF for point at (7,8,8); SMR recon- struction. Z-SECTION 8: X: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Y/ O 0 0 0 O O 0 0 O 0 0 0 0 0 O O 0 1 O 0 0 0 O 0 0 0 O O O 0 0 0 0 0 2 0 O O 0 0 O 0 0 0 0 0 0 0 O 0 O 3 0 O O 0 0 0 O *0 0 0 O 0 0 O 0 0 4 O 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 O 0 0 0 0 0 0 0 O 0 O O 0 0 0 6 O 0 0 0 O 0 0 0 0 0 0 O O O 0 O 7 0 0 O O 0 0 0 0 O 0 0 O O O 0 O 8 0 0 O 0 0 0 10 30 0 ‘0 0 O 0 O 0 0 9 O 0 0 O 0 0 O 10 0 O 0 0 O O 0 0 10 O 0 O 0 0 0 0 0 0 O 0 0 0 O 0 O 11 0 O 0 O 0 O 0 0 O 0 O 0 0 0 O 0 12 O O 0 O 0 0 0 0 0 0 0 0 0' 0 O O 13 0 0 O O O 0 0 0 0 0 0 O O 0 0 O 14 O O 0 0 0 0 0 0 0 O 0 0 0 0 0 0 15 0 O O 0 0 O 0 O 0 0 0 0 O O O O 42 (Y = 12 through 15 all zeroes) Poor—longitudinal-resolution PSF for point at (7,8,8); SMR reconstruction. Table 4-6. Z-SECTION 8: 9 10 11 12 13 14 15 8 2 1 0 X 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 00000000101000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 0123145678901 11 9 10 ll 12 13 14 15 8 l 2 0 Z-SECTION 9: 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 oooooooomooo 000000000000 oooooooomooo 000000000000 000000000000 000000000000 000000000000 000000000000 0112345678901 11 9 10 11 12 13 14 15 8 2 1 O Z-SECTION 10: x. 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 oooooooomooo 000000000000 000000000000 oooooooomooo 000000000000 000000000000 000000000000 000000000000 4000000000000 0112345678901 11 43 constant-valued as is that of the MRR, but has a maximum at the true location of the ideal reflector. This suggests that for the summation rule some form of post-processing of the reconstructed image may be possible in order to minimize the image-degrading effects of the PSF. The PSF for the sum of squares rule (SSR) under resolution- matched conditions is shown in Table 4-7. The PSF under conditions of poor longitudinal resolution appears in Table 4-8. The actual loca- tion of the point reflector in each case is at (7,8,8); the point's echo has a measured value of 10 from each View. The character of the SSR PSF is similar to that of the summation rule PSF. Again, the ability to compensate for the PSF in some way is suggested by the fact that the PSF has a maximum at the true location of the ideal point reflector. Table 4-7. Resolution-matched PSF for point at (7,8,8); SSR recon- struction. Z-SECTION 8: X: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Y/ O 0 0 0 0 O O 0 O 0 0 O 0 0 O 0 0 1 0 0 O O O O 0 0 O O O O O O 0 0 2 0 0 0 O O O O 0 0 0 0 0 0 0 0 O 3 0 0 0 0 0 0 0 0 0 0 O 0 0 0 0 O 4 0 O 0 0 O O 0 0 0 0 0 0 O 0 0 0 5 O 0 0 0 0 0 0 O 0 O 0 0 O 0 0 0 6 0 O O 0 0 O 0 0 O 0 O 0 0 O 0 0 7 0 O 0 0 0 0 0 0 0 0 0 0 0 O O 0 8 0 0 0 0 O O 100 300 0 O 0 0 0 0 O 0 9 0 0 O 0 0 0 0 100 0 0 0 0 0 0 0 0 10 O 0 0 0 0 O O O 0 0 0 0 0 O 0 0 11 0 0 0 0 0 0 0 0 0 0 0 0 0 O 0 0 12 0 O 0 0 0 0 0 O 0 0 0 0 0 0 O O 13 0 0 0 0 O O 0 O O 0 0 0 0 0 O O 14 0 0 0 0 0 0 0 0 0 0 0 O O 0 0 O 15 0 0 O 0 0 O 0 0 O 0 0 0 0 0 0 0 44 Poor-longitudinal-resolution PSF for point at (7,8,8); Table 4-8. (Y = 12 through 15 all zeroes) SSR reconstruction. Z-SECTION 8: 9 10 11 12 13 14 15 8 2 1 0 X II V. 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 00 31 oooooooomooo 1 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 012345678901 .11 9 10 11 12 13 14 15 8 2 1 0 Z-SECTION 9: X II V. 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 oooooooomooo .I. 000000000000 00 0 11... .l. 000000000000 OOOOOOOOWOOO 1... 000000000000 000000000000 000000000000 000000000000 000000000000 012345678901 11 Z-SECTION 10: 9 10 11 12 13 14 15 8 2 1 0 .II v. 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 OOOOOOOOWOOO 1 000000000000 000000000000 0 0 0 .l. 1 1|. 000000000000 oooooooomooo .I. 000000000000 000000000000 000000000000 000000000000 000000000000 012345678901 11 5. PERFORMANCE OF RECONSTRUCTION METHOD The performance of a reconstruction method may be evaluated empi- rically through the use of data measured on a "phantom". A phantom is an actual object whose physical dimensions and characteristics are known to the experimenter. The goodness of a reconstruction technique may be judged in pr0portion to its ability to reconstruct an image which reliably manifests the properties of the phantom which the tech- nique is designed to image. It is important to take into account the limitations of the measurement or data acquisition apparatus in the evaluation. The ELECTRA program was tested on data obtained from a simple phantom. A description of the phantom and data acquisition system, and a discussion of the test results, follow. 5.1 Test data. The phantom from which the test data were measured consists of two square plexiglas sheets, roughly six inches on a side, held paral- lel to each other about two and one-half inches apart by four bolts near the corners of the sheets. In addition, plastic monofilament fishing line is strung through holes in the sheets across the space between the sheets and perpendicular to the faces of the sheets. 45 46 The measurement apparatus consists of a 3 MHz nominal-center- frequency transducer, with face diameter 10 mm, which rides on a screw rod over the surface of a water bath (in which the object of interest lies -- the water is the propagating medium for the ultrasound). The face of the transducer is submerged in the bath, and a step motor on the screw rod allows the transducer to be scanned across the bath in discrete user-programmed steps. After the transducer is moved to a new position along the rod, it is excited to produce the interrogating pulse, then switched to the listening mode to record echoes. The echoes are full-wave rectified and digitized to 8 bits at 5 MHz for a total of 2048 samples. The digitized data are sent to a Varian V75 minicomputer for storage on magnetic tape. The V75 controls the acti- vity of the entire measurement apparatus. For the purpose of generating a data set with which to test the ELECTRA software, echo data were used from a scan performed on the phantom by the measurement device in which the transducer was scanned down a line running parallel to and between the plexiglas sheets. In this way, the bolts and fish lines traversing the space between the sheets were scanned in cross-section. The step size used during the scan was 3 mm. The "longitudinal" data values for each scan position were integrated such that the longitudinal sample extent was similar to the 3 mm lateral sample extent. The echo samples are digitized to values between zero and 255, but in order to keep the numbers small in this test reconstruction these values were scaled to a range of zero to 10. The reconstruction region of interest was selected to contain one bolt and one fish line. The bolts are 3.5 mm in diameter, representing 47 large non-Rayleigh scatterers; the fish line is 0.5 mm in diameter, approximating a Rayleigh scatterer for 3 MHz ultrasound. Because the bolts and fish line are cylindrical, their cross-sectional scattering is isotropic; any scan made across a cylindrical object perpendicular to its axis may be assumed to be equivalent to any other such scan of the same object. Utilizing this assumption, the data obtained from the vertical view of the bolt and fish line phantom were used to simu- late views of the phantom from other angles, by making appropriate changes to the data to account for the different relative positions of objects due to changes in their angle of projection (view). The data set resulting from these manipulations simulates data obtained from scanning a region of interest through which two "anoma- lies" run parallel to each other. A schematic view of the problem is shown in Figure 5-1. This region is assumed to have been scanned from three views: vertical -- in which the scan lines run perpendicular to the axes of the bolt and fish line; left and right -— both at 45 degree angles, in which the scan lines run parallel to the axes of the bolt and fish line. In the vertical view, the data values for the i-th element in all scan lines are equal. In the right and left views, all data values in each scan line are equal. The reconstruction region is 16 x 16 x 16 voxels; the scans are 16 x 16 voxels. Table 5-1 shows the data obtained in the vertical view. Note that this figure and all subsequent figures in this section show "y-sections" in which the bolt and fish line are seen in cross—section. All y-sections are equal and any one y-section demonstrates the pro- perties of the reconstruction technique in any dimension. 48 left vertical right iew View vie x -————-'. observation A height 5 ' 1’ " units /// / I I / I /' bolt | I’ z I / ’/ fish line / Figure 5-1. Schematic drawing of test reconstruction problem. 49 Two things should be noted in Table 5-1: (1) The bolt (located at (x,y,z) = (4,y,2)) casts a long "shadow" of echoes. Some echoes come from reflection at the "front" edge of the water-bolt interface and some from the "back" edge (bolt-water). At (4,y,9) a second large echo is measured. The transducer is 5 units above the top (2 = 0) of the R8 in this problem, that is 20 = -5. The bolt is at z = 2, 7 units from the transducer face. The second bolt echo at (4,y,9) is 14 units from the transducer face. This echo is due to ultrasound following the path: transducer (origin) to bolt (reflection) to transducer (reflection) to bolt (reflection) to trans- ducer (measurement). This multiple reflection path represents a com- mon source of artifacts in ultrasonic imaging [16]. (2) The images of the objects show considerable lateral spread. This is because the bolt, measuring 3.5 mm in diameter, and the line, Table 5-1. Vertical view test data. Y-SECTION 1: x: o 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 z/ o o o o o 0 o o o o o o o o o o 0 1 o o 0 o o 1 o 0 o o 0 o o 0 o o 2 o 1 6 10 10 9 3 1 o o o o o o o 0 3 o 2 1 9 10 6 2 1 o o o o 0 o 0 o 4 1 1 1 6 9 2 o 1 o 0 0 o o o o o 5 o o 1 6 9 1 o o o o o o o o o o 6 o o o 4 8 1 o o o o 0 o 0 o o o 7 o 0 o 2 6 1 o o o o 0 o o o o o 8 o o o 1 2 1 o o o o 1 1 1 o o 0 9 0 o 1 7 9 2 o 0 o 1 5 7 2 1 o o 10 0 o 0 6 9 1 o 0 o o 1 o o 0 o 0 11 0 o 0 3 8 1 0 0 0 0 o o o 0 o o 12 o o o 4 4 o o o o o o o o 0 o o 13 o o o 2 3 o o 0 o o , 0 o o o o o 14 0 o o 1 1 o o o 0 o o o o o o o 15 o o o o 1 o o o 0 0 0 o o o 0 0 50 measuring 0.5 mm in diameter, remain in the transducer beam for sever- al measurements. The transducer is moved along in 3 mm steps, but its face (and its main, near-field beam) measures 10 mm in diameter. In addition, any unfocussed ultrasound transducer similar to the one used to obtain these data has "side lobes" to its beam which extend beyond the cylinder defined by the face of the transducer. Energy re- flected from the side lobes would be interpreted as detection of a re- flector in the main beam. Table 5-2 shows the voxels for which non-zero data were obtained in the vertical view. This figure is produced by binary windowing of the data values such that all values greater than or equal to a threshold (1.0 in this case) are displayed as "@" and all values less than the threshold are displayed as Table 5-2. Binary—windowed vertical view test data; threshold = 1. Y¥SECTION l: X: 0123456789ABCDEF Z ......@....... .... .@@@@@@@........ .@@@@@@@........ @@@@@@.@... ..... ..@@@@.......... ...@@@.......... ...@@@.......... ...@@@....@@@... ..@@@@...@@@@@.. ...@@@....@..... ...@@@.......... ...@@........... ...@@........... ...@@..... ..... . ....@............. WNUOW>OmN®MbWNHO\ 51 The data obtained for the left view are shown in Tables 5-3 and 5-4. Note that from this view the bolt shadows the fish line. Tables 5-5 and 5-6 show the right view data. Values of "-1" in these tables indicate voxels for which no measured data were obtained. 5.2 Maximum-replacement rule reconstruction. Table 5-7 shows the reconstruction obtained under the maximum- replacement voxel value modification rule (MRR). The bolt is at (4,y,2) and the fish line is at (11,y,9). Binary windowing of the reconstruction at varying thresholds manifests the character of the MRR results. Selecting an arbitrary threshold of 5.0, one sees in Table 5-8 that the signal-to-noise ratio is extremely poor for the MRR. Increasing the threshold to 7.0 allows one to isolate the fish line (see Table 5-9), but the true shape and size of the bolt is not yet clear. In order to isolate the bolt, one must increase the threshold to the maximum level of 10.0 (no larger data value will appear under the MRR for this test case). This is shown in Table 5-10. The bolt is isolated at the expense of "losing sight" of the fish line. While gray-scale imaging would allow an observer to see the bolt and the line, it would also allow the observation of the strong secon- dary bolt reflection. Binary windowing of the reconstruction is use- ful for illustrating a problem common in ultrasonic imaging: An arti- fact of a strong reflector (the secondary bolt echo) may be given greater reflector-presence credibility than the direct echo of a weak reflector (the fish line). This observation suggests a rule-of-thumb 52 Left view test data. Y-SECTION 1: 9 10 11 12 13 14 15 8 2 1 0 X / Y -1 -1 -1 4000000000004 -1 -1 -1 0 -1 0000000000000 -1 -1 0 -1 0 -1 0 0 0 -l 0 1 -l 1.0.4.000 -1 0 0 0 -1 0000001012000000 0101462100000000 11146914004001.1000 00.192100000000011. 0.1234567890112345 111.111... Binary-windowed left view test data; threshold = 1. Table 5-4. X:0123456789ABCDEF Y-SECTION 1: z/ .@...1. ..@@@@. ...@@@.@.... 0 1 2 3 O... .@@... ...@@@@@@.... ..@@.@@.@. ...@.@@@@@... .......@@.@.. A.§.Av7.R.Q,Aununununuw. 53 Right view test data. Table 5-5. Y-SECTION l: 9 10 11 12 13 14 15 8 2 1 0 X I Y —1 -1 -1 -1 -1 -1 -l -1 -1 0 -1 0 -l —l -1 -1 1111111110001111 _.___.._. .... 1.111111001100011 ...—... ._ 1111111012100101 ______. . . _ 57000000 -1 -1 -1 —1 -1 0 0 1 _ . 0000000000000 -1 -1 0 -1 -1 0 0 0 0013100000000000 -1 -1 9 9 1 -1 0 -1 -1 6 -1 9 1 -1 61100000000 -1 -1 -1 6 8 -l -l -1 Jul—12200000000 01214456789019.3145 111111 Binary-windowed right view test data; threshold = 1. Table 5-6. F E D C B A 9 8 7 0.6 15 I4 ms 12 m1 Em 48X Y Z/ 0 @.. .@@@@@.... .@.@@.@......... 4 @@@@@@@......... 5 @@@..@. 1 2 3 6 @@.... ...@@. @ @ 7 .@@@.. ....@. ...@@... 89ABCDEF 54 Maximum replacement rule test reconstruction. Table 5-7. Y-SECTION 1: 9 10 11 12 13 14 15 8 0000000000000000 0000000000000000 0000000001100000 0000000012100000 0000000157000000 0000001115200000 0000010201000000 0.101462100000000 0.216181000000000 0.136910000000000 1.090211112110000 1 1 0 09986299814311 1 0 10 0 96642176314210 0 1 0 1 0061910001000000 0016861100000000 0000212200000000 val/nu119.2.4.RJAV7.Ruo,nv1.9.14A.§. 11111111 2 Binary-windowed MRR test reconstruction; threshold = 5. Table 5-8. P. p“ D no R. _A o, no 7 .06 1.RJ .4 ma T.9. m1 Em 4vA Y ...@@.... .@@... ... .@....... .@@....... ....@@... ooooooo@oooo .@@@@........ .@.@@@@@...... .@@@@.@. .@.@@. ....@. ...@. ...@@.. [Inu119.1.h.fiafiu7.RVOIAuRapunuwuwr 7g 55 Binary—windowed MRR test reconstruction; threshold = 7. Y-SECTION 1: X: z/ WHUOW>O®N¢UIJ>WNHO 0123456789ABCDEF .. .@@........:. ...@@@... ....... .............@@@ .@@.@.@... ...... ....@..@.. ......@ ...@@......@...: Binary-windowed MRR test reconstruction; threshold = 10. Y-SECTION 1: X:0123456789ABCDEF z/ Hinic3633!>~¢>03\40\U1¢~u>n>hic> 56 for the interpretation of ultrasonic images: Consider the maximum- valued reflections first and assume that they represent actual reflec- tors. Next consider lesser-valued reflections in order (proportional to their value). If a weaker reflection can be eXplained with respect to existing "actual" reflectors and the observer's knowledge of the reconstruction domain, assume that they are not the product of real reflectors. If weak reflections cannot be explained, then assume that they are the product of a weak reflector. 5.3 Summation rule reconstruction. The reconstruction obtained under the summation rule (SMR) is shown in Table 5-11. Notice that the location of the bolt (4,y,2) and Table 5-11. Summation rule test reconstruction. Y-SECTION l: x: o 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 z/ o o o 0 o o o 1 o o 0 o o o o o o 1 0 o o 0 10 17 1 2 2 0 o o o o o 0 2 o 1 12 20 4o 27 4 2 o o o o o o o o 3 o 8 1 3o 22 25 11 7 1 o o o o o o 0 4 3 13 11 9 11 6 11 1 4 0 o o o o 0 0 5 1 6 2 8 9 3 1 8 8 1 o o o o o 0 6 3 1 0 4 8 1 o 1 2 0 1 0 o 0 0 o 7 3 1 o 2 6 1 o o 1 2 1 1 o o o o 8 o 0 o 1 2 1 o o o o 3 6 3 o 0 0 9 o o 1 7 9 2 o o o 1 5 19 4 2 o o 10 o o 0 6 9 1 o o o o 3 o 1 1 o o 11 o o 0 3 8 1 o 0 o o 0 o o 0 0 o 12 0 o o 4 4 o 0 0 o o 0 o o 0 o 0 13 o 0 o 2 3 o 0 0 0 o 0 0 0 o o 0 14 o o o 1 1 o o o o 0 o o o o o 0 15 o 0 o o 1 o o 0 o o 0 o o o o o 57 the fish line (11,y,9) are correctly indicated by the local maxima of the reconstruction array. Notice also that the secondary bolt echo artifact has a value less than that of the fish line. This second observation is illustrated in Table 5-12, in which a binary window threshold of 15.0 is used. This threshold is the SMR equivalent of the MRR threshold of 5.0 (3 views x 5 = 15); thus, Tables 5-12 and 5-8 may be compared. The SMR is clearly a rule of superior properties to the MRR. By increasing the window threshold to 19.0, as shown in Table 5-13, one may achieve an acceptable representation of the structure in the RS. The fish line is correctly located and is one voxel in cross- sectional dimension. The bolt is correctly located and its cross- sectional extent is as well defined as may be expected under the limi- tations of the measurement and reconstruction method: the bolt is 3.5 mm is diameter and therefore may be expected to span from 2 to 4 of the 3 mm-on-a-side voxels; the reconstruction error of i1 voxel along any axis can account for the 2 x 3 voxel cross-section recon- structed here. 5.4 Sum of squares reconstruction. Table 5-14 shows the reconstruction obtained under the sum of squares rule (SSR). Like the SMR reconstruction, the SSR reconstruc- tion correctly locates the bolt and fish line by its local maxima. Binary windowing with a threshold of 75.0 (equivalent to the SMR threshold of 15.0 and the MRR threshold of 5.0 -- 3 views x 52 = 75) indicates that while large signals are amplified in proportion to 58 Table 5-12. Binary-windowed SMR test reconstruction; threshold = 15. Y-SECTION 1: X:0123456789ABCDEF z/ MMUOUU>KOCDNOLnwaHO Table 5-13. Binary-windowed SMR test reconstruction; threshold = 19. Y-SECTION 1: x:0123456789ABCDEF Y/ .3.@@@.......... .............@@@ WMUOWIDKOWNO‘M-bWNF-‘O 59 their value by the SSR, so are large noise values (see Table 5-15). The secondary bolt echo artifact shows up at (4,y,9) and (4,y,A). Increasing the threshold to 115, as shown in Table 5-16, allows one to isolate the fish line and greatly reduce the noise surrounding the bolt. A reconstruction of the bolt similar to that for the SMR with threshold 19 may be achieved at a SSR threshold of 170 (see Table 5-17), but at the cost of losing the fish line. 5.5 General comments on performance. The summation (SMR) and sum of squares (SSR) rules have a dis- tinct signal-to-noise advantage over the maximum-replacement rule (MRR). This is due to the reinforcement of multiply measured, locally Table 5-14. Sum of squares rule test reconstruction. Y-SECTION 1: X: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Y/ 0 O 0 0 0 O 0 1 0 0 0 0 O 0 0 0 0 1 0 0 O 100 136 1 4 2 0 0 0 0 O O O 0 2 0 1 72 181 400 242 10 2 O 0 0 0 O O 0 0 3 0 4O 1 251 171 216 49 37 1 O 0 0 O 0 O 0 4 5 81 82 41 84 12 82 l 16 0 0 0 0 0 0 0 5 1 36 2 38 80 3 1 64 4O 1 0 0 0 O O O 6 5 1 0 16 64 1 O 1 4 0 l 0 0 0 0 0 7 5 1 0 4 36 l 0 0 1 4 1 1 0 0 0 0 8 O 0 O 1 4 1 O 0 O 0 3 26 3 0 O 0 9 0 0 l 49 8O 4 0 0 0 l 25 115 8 2 O 0 10 O 0 0 36 80 1 0 0 O 0 5 0 1 1 O O 11 0 0 O 9 64 l 0 O 0 O O O O 0 0 0 12 0 O 0 16 16 0 0 O O O O O 0 0 0 0 13 O 0 0 4 9 0 0 0 0 0 0 0 O 0 O 0 14 0 O O 1 1 O O 0 0 O 0 O 0 O 0 0 15 O 0 0 0 1 0 0 O O 0 0 0 0 O O 0 60 Table 5-15. Binary-windowed SSR test reconstruction; threshold = 75. Y-SECTION 1: :0123456789ABCDEF Z m131c1c>u:>-m>a>xloxuan~u>n>hac>\e>< Table 5-16. Binary-windowed SSR test reconstruction; threshold = 115. Y-SECTION 1: :0123456789ABCDEF Z . ....... ..@.. m1u1c1c:u:m-u>a:~JO\U1¢~uJhawuca->4 61 co-occuring echo values by their sum. Weak anisotropic reflectors would not benefit from this reinforcement, but these reflectors would be difficult to image under any of the rules discussed in this thesis. The SSR appears to suffer from its measured-value weighting scheme (ie. squaring the value before summing it) in that it is possi- ble for single large non-zero measurements of artifacts to attain a greater value in their square than the sum of squares of multiple small non-zero measurements of weak, but real, reflectors. For this reason, and in light of the results of the test reconstruction, one may wish to declare the SMR to be the best reconstruction rule among those discussed. However, such a decision may be premature, based as it is on limited evidence. Table 5-17. Binary-windowed SSR test reconstruction; threshold = 170. Y-SECTION 1: X:0123456789ABCDEF z/ H1n1Cir:u1>~¢>a>\40\uwc~uan>-c> 6. SUMMARY AND CONCLUSIONS A method has been described for the reconstruction of three- dimensional images from two-dimensional ultrasonic data. The method may be implemented in modular computer software, making it easy to adapt to changes in measurement and reconstruction geometry, and to varying post-processing requirements. In additibn, the method pro- vides for incorporation of newly-acquired data into the existing re— construction as soon as the data become available. Thus, the recon- struction may be viewed as it evolves. The reconstruction method has two steps. First, the spatial coordinates for a measured data point are transformed from the mea- surement frame of reference to the reconstruction frame of reference. Second, the newly measured data value is combined in some way with the existing value at the reconstruction-frame location. The point spread function (PSF) of this reconstruction method was investigated empirically. When the spatial resolution of the measure— ment apparatus is equal to or better than that of the reconstruction array, the spreading of an image of a point reflector is minimal and is due to quantization error produced by rounding off computed recon- struction space coordinates to integer values used for addressing ele- ments in the reconstruction array. When the spatial resolution of the measurement apparatus is not as good as that of the reconstruction array, the PSF becomes more extended and takes on a complex shape. 62 63 The nature of point image degradation described by the PSF is dependent on the technique used to modify existing reconstruction array values by combination with newly measured data. Three tech- niques were investigated: replacement of the existing value with the maximum of the existing and newly measured values; accumulation of the sum of measured values; accumulation of the sum of squares of measured values. The maximum replacement rule has the poorest performance due to the property of granting equal reflector presence credibility to all points for which a given echo amplitude has been recorded at least once. The summation rule has the advantage of reinforcing data values for repeated (confirming) measurements of a reflection at a given location. The sum of squares rule shares that advantage and further reinforces large measured values by weighting them by them- selves. The superior performance of a summation or integrate rule with respect to that of a maximum replacement or peak rule has been con- firmed by concurrent, independent research [17]. In most digital pro- cessors, summation (and, often, summing squares) may be performed in less time than required for the computation of a maximum. All three rules require only as much memory as is needed for the reconstruction array, and allow the user to view the reconstruction as it evolves. Post-processing of the reconstruction array may be used to en- hance the information portrayed in the image. Binary windowing con- sists of selecting a threshold value above which one considers data to represent signal and below which one considers data to represent noise. Binary windowing was used in this study to characterize reconstruction performance. Human viewers of ultrasonic images perform visual 64 step-wise binary windowing, interpreting strong signals first, then working toward the interpretation of weaker signals. Machine post—processing by surface-seeking algorithms has been demonstrated to be useful in clinical applications [18—20]. Ultra- sonic scattering is caused by changes in the physical characteristics of the propagating medium. Since these changes are typically indica- tive of material interfaces, discovery and display of surfaces may be an effective aid to the interpretation of a reconstruction. Ultrasonic measurement is susceptible to artifacts resulting from the effects of scattering phenomena (diffraction, refraction, secon- dary reflection). Strong reflection signals may be observed for loca- tions at which no real reflector exists. No amount of reconstruction post-processing can be as effective in aiding image interpretation as can be knowledge of the general structure of the region of interest. In the materials testing context, knowledge of the location and shape of "intentional" structures such as drill holes or fasteners can pro- vide the image viewer with landmarks by which to judge the validity (location and echo intensity) of unexpected reflections possible due to defects. The evaluation of the reconstruction method presented here has been based on data obtained from a phantom of simple structure under somewhat idealized conditions. Future work in this area should in- clude evaluation of the method's performance using data obtained from a test object in an actual materials testing application. In particu— lar, further investigation is warranted into the effect on reconstruc- tion performance of the number and angle of oblique views, and the presence or absence of structural symmetry and anisotrOpic reflectors 65 within the volume of interest. With respect to medical applications of this work, the recon- struction method may not be applicable generally to three-dimensional diagnostic imaging. This is due to the problem of object-of—interest motion in clinical imaging and the subsequent requirement of a very rapid data collection scheme. That condition is not likely to be satisfied by the method for obtaining the two-dimensional data requir- ed by this reconstruction technique. However, it has been demonstra- ted in this thesis that the method of modifying a reconstruction voxel value by the sum of measured values has associated with it a point spread function more amenable to image enhancement by post-processing than does the method of modification by replacement with the maximum observed value. Effort should be made to evaluate the performance of integrate or summation modification rules in clinical imaging systems, which commonly use a maximum-replacement or peak rule due to its con— venient implementation in analog hardware [17]. Further development of digital techniques and their application to ultrasonic imaging will continue to enhance the utility of this nondestructive, noninvasive imaging method. APPENDIX APPENDIX ELECTRA PROGRAM LISTING ********************************************************************** TITLE: ELECTRA AUTHOR: D.A. GIFT DATE: 05 MAY 1980 PURPOSE: PRODUCE 3-D RECONSTRUCTIONS FROM 2-D BACKSCA'I'I‘ERED ULTRASOUND DATA. ******* ********************************************************************** DESCRIPTION OF VARIABLES ANGLE - INPUT PARAMETER FOR EACH SLAB; ANGLE OF TRANSDUCER AIM WITH RESPECT TO zeAXIs 0F RECONSTRUCTION SPACE. ATSIGN - "@"; BINARY-WINDOWED OUTPUT SYMBOL. BINWIN - 16 x 16 x 16 ARRAY USED FOR OUTPUT OF RECONSTRUCTION ARRAY IN VARIOUS WAYS WITHOUT ALTERING CONTENTS OF RECONSTRUCTION ARRAY ITSELF; MUST BE DIMENSIONALLY EQUAL TO "RMATRX"O DATA - 256-ELEMENT ARRAY USED As INPUT BUFFER FOR MEASURED VOXEL DATA (THIS PROGRAM ASSUMES A 16—LINE BY 16 RASTER UNIT SCAN FOR A TOTAL OF 256 SAMPLES/SLAB. DAUG - 3 x 4 MATRIX FOR AUGMENTED MATRIx.OF DIRECTION COSINES USED IN COORDINATE TRANSFORMATION. DEPTH - INPUT PARAMETER FOR EACH SLAB; DEPTH OF SLAB FROM TRANSDUCER. DOT - "."; BINARY-WINDOWED OUTPUT SYMBOL. DSPLAY - USER-INPUT SELECTION OF RECONSTRUCTION- SPACE AXIS FOR OUTPUT SECTIONS. FCBB - 13-ELEMENT ARRAY; FILE CONTROL BLOCK FOR DISC FILE FROM WHICH INPUT MEASURED DATA ARE READ (FILE NAME - "URDATA"); PECULIAR TO LOCAL I/O CONTROL UNDER VARIAN VORTEX II OPERATING SYSTEM. IRDIM - LENGTH OF EACH SIDE OF THE RECONSTRUCTION ARRAY; HERE SET EQUAL TO 16. IX - RECONSTRUCTION-SPACE XHAXIS COORDINATE. IY - RECONSTRUCTION-SPACE YHAXIS COORDINATE. IZ - RECONSTRUCTION-SPACE ZeAXIS COORDINATE. KSLAB - COUNTER FOR NUMBER OF SLABS INPUT. LABEL - 16-ELEMENT ARRAY OF SYMBOLS FOR OUTPUT TABLES. a-a-aa-a-a-a-a-a-a-a-a-a-a-a-aa-a-anx-a-a-a-a-W-x-a-a-a-a-a-a-a-a-a- 66 i * i i * * * fl * fl * fi * i * i i i i i t i * i i i * i i i % & t i i i i * MPCODE NDATA OUTLUN RMATRX RMMRRT RMSMRT RMSSRT THRESH VOXNUM X0 XSECTN YO YSECTN ZO ZM ZSECTN 67 INPUT PARAMETER FOR EACH SLAB; MEASUREMENT PLACEMENT OR VIEW CODE: V=VERTICAL, F=FRONT, B=BACK, R=RIGHT, L=LEFT. NUMBER OF MEASURED VOXELs/SLAB; SET EQUAL TO 256 HERE; EQUAL TO DIMENSIONALITY OF "DATA" ARRAY. COMPUTER SYSTEM LOGICAL UNIT NUMBER OF OUTPUT DEVICE; THIS PROGRAM ASSUMES A TEKTRONIX 4014-1 GRAPHICS TERMINAL. 16 x 16 x 16 ARRAY; STORAGE FOR RECONSTRUCTION ARRAY. CONTAINS MNEMONIC FOR INTERACTIVE SELECTION OF MAXIMUM-REPLACEMENT RULE. CONTAINS MNEMONIC FOR INTERACTIVE SELECTION OF SUMMATION RULE. CONTAINS MNEMONIC FOR INTERACTIVE SELECTION OF SUM OF SQUARES RULE. USED IN INTERACTIVE OUTPUT SELECTION. INPUT PARAMETER FOR EACH SLAB; SLAB IDENTIFICATION NUMBER; BOOKKEEPING PURPOSES ONLY. INTERACTIVE THRESHOLD VALUE FOR BINARY-WINDOWED OUTPUT. NUMBER OF CURRENT VOXEL IN CURRENT SLAB. RECONSTRUCTION-SPACE x.COORDINATE OF CENTER OF SCANNING APERTURE. MEASUREMENT-SPACE x COORDINATE. CONTAINS MNEMONIC FOR INTERACTIVE SELECTION OF x-SECTION OUTPUT TABLES. RECONSTRUCTION-SPACE Y COORDINATE OF CENTER OF SCANNING APERTURE. MEASUREMENT-SPACE Y COORDINATE. CONTAINS MNEMONIC FOR INTERACTIVE SELECTION OF Y-SECTION OUTPUT TABLES. RECONSTRUCTION-SPACE z COORDINATE OF CENTER OF SCANNING APERTURE. MEASUREMENT-SPACE z COORDINATE. CONTAINS MNEMONIC FOR INTERACTIVE SELECTION OF Z-SECTION OUTPUT TABLES. ********************************************************************** REAL DATA(256), DAUG(3,4) INTEGER SLAB, VOXNUM, LABEL(16), OUTLUN, DOT, ATSIGN INTEGER BINWIN(16,16,16), XSECTN, YSECTN, ZSECTN, S, DSPLAY INTEGER FCB8(13)/2*O,' ',4*0,'UR','DA','TA’,3*0/ COMMON /RECON/ RMATRX(16,16,16), IRDIM, KSLAB DATA NDATA/256l, RMMRRT/’MRRT'l, RMSMRT/‘SMRT'I DATA RMSSRTI'SSRT'l, N/’N'/, DOT/'.’/, ATSIGN/'@’/ DATA XSECTN/'x’/, YSECTNI’Y’I, ZSECTN/'Z'l, s/'s’/ DATA (LABEL(I),I=1,16)/'O’,‘1',’2',’3',‘4','5','6','7',’8', + ’9','A”'B’"C"'D”’E”’F’/ 68 * ********************************************************************** i: 1000 FORMAT(’O’,A1,’-SECTION', I3, ':’) 1010 FORMAT(1X,I2,ZOI3) 1020 FORMAT(‘1') 1030 FORMAT(//,’OEXECUTION TIMEa',F8.O,’ MILLISECONDS,',/, + ' LESS',F8.O,' MILLISECONDS (ESTIMATED) DISC ACCESS TIME,',/, + 'OYIELDS EFFECTIVE TIME FOR RECONSTRUCTION -',F8.0, + ’ MILLISECONDS.’) 1040 FORMAT(zx,A1,’:',Iz,15I3,/,1x,A1,'/') 1050 FORMAT(A4,1x,Iz) 1060 FORMAT('OINPUT CHOICE OF RECONSTRUCTION TECHNIQUE:’,/, T5,' FOR MAXIMUM REPLACEMENT RULE;',/, T5,' FOR SUMMATION METHOD;’,/, T5,' FOR SUM OF SQUARES METHOD.',/, 1X,'AND OUTPUT UNIT # -- ,.') 1070 FORMAT(’ MAXIMUM REPLACEMENT'RULE USED FOR RECONSTRUCTION.’) 1080 FORMAT(’ SUMMATION METHOD USED FOR RECONSTRUCTION.’) 1090 FORMAT(' SUM OF SQUARES METHOD USED FOR RECONSTRUCTION.’) 1100 FORMAI(’0DO YOU WANT BINARY WINDOW OUTPUT...(Y OR N)’) 1110 FORMAT(A1) 1120 FORMAT(’OINPUT THRESHOLD VALUE...(F10.O)') 1130 FORMAT(F10.0) 1140 FORMAT(IX,A1,’-SECTION ',A1,T22,A1,'-SECTION ’,A1,T42,A1, + '-SECTION ',A1,/, + 2x,3(16A1,4X)) 1150 FORMAT(lx,3(17A1,3X)) 1160 FORMAT(lx,A1,’-SECTION ',A1,9x,’MOD. RULE: ’,A4, + ' THRESHOLD -’,F5.0,/,2X,16A1) 1170 FORMAT(lx,17A1) 1180 FORMAT(l/l) 1190 FORMAT(l) 1200 FORMAT(//) 1210 FORMAT('OINPUT RECONSTRUCTION DISPLAY SECTION AXIS: (A1)’, + /,' ( - SKIP DISPLAYY) + ... + .9. a: ********************************************************************** CALL V$0PEN(8 , 98,FCBB, O) WRITE (60, 1060) READ (60, 1050) RMETHD,OUI‘LUN CALL TIME (MINS , MILLIS) IRDIM-1 6 KSLAB-O * INITIALIZE RMATRX T0 "NO DATA". *- DO 100 I=1,IRDIM DO 100 J-l,IRDIM DO 100 R-1,IRDIM RMATRX(I,J,K)--l. 100 a! 69 CONTINUE ********************************************************************** a-ara-a-x- READ SLAB DATA FROM MASS STORAGE; GENERATE MS VOXEL COORDINATES; DRIVE RECONSTRUCTION. ********************************************************************** * * READ SLAB NUMBER, DEPTH FROM TRANSDUCER, TRANSDUCER POSI- * TION IN RS, TRANSDUCER MEASUREMENT PLACEMENT CODE, ANGLE OF * TRANSDUCER AIM (RADIANS), AND MS VOXEL DATA. * 10 READ(8,END=6OO) SLAB, DEPTH, x0, YO, zo, MPCODE, ANGLE, + (DATA(R),R-1,NDATA) KSLAB=KSLAB+1. *- * FORM AUGMENTED DIRECTION-COSINE MATRIX. *_ CALL DIRCOS(MPCODE,ANGLE,DAUG) DAUG(1,4)'XO DAUG(2,4)-YO DAUG(3,4)-20 *— * GENERATE VOXELS BELONGING TO CURRENT SLAB; * CARRY OUT RECONSTRUCTION. *- DO 500 VOXNUM=1,NDATA CALL VOXGENKDEPTH,VOXNUM,XM,YM,ZM) * * APPLY SPATIAL TRANSFORM TO CURRENT MS VOXEL. *= CALL TRNSFM(XM,YM,ZM,DAUG) *_ * CALL RECONSTRUCTION ROUTINE WITH CURRENT VDXEL. 1*, IZ'INT(ZM+.5)+1. IF(IX.LE.O .OR. IY.LE.O .OR. IZ.LE.O) GOTO 500 IF(IX.GT.IRDIM .OR. IY.GT.IRDIM .OR. IZ.GT.IRDIM) GOTO 500 IF(RMETHD.EQ.RMMRRT) CALL MRRT(IX,IY,IZ,DATA(VOXNUM)) IF(RMETHD.EQ.RMSMRT) CALL SMRT(IX,IY,IZ,DATA(VOXNUM)) IF(RMETHD.EQ.RMSSRT) CALL SSRT(IX,IY,IZ,DATA(VOXNUM)) 500 CONTINUE NEXT SLAB. GOTO 10 7O * IF NO MORE SLABS. . . 1.- 600 CONTINUE CALL T IME (MINE , MILLIE) RT IME- ( (MINE-MINS )*1 . 2E4+(MILLIE-MILLIS ) )*5 . DTIME=KSLAB*55. ET IME=RT IME-DT IME CALL V$CLOS (8.0) * PRINT RECONSTRUCTION MATRIX. *_ 601 WRITE(6O,1210) READ(6O,1110) DSPLAY IF(DSPLAY.EQ.S) GOTO 700 D0 603 I-1,IRDIM DO 603 J-1,IRDIM DO 603 K-1,IRDIM IF(DSPLAY.EQ.ZSECTN) BINWIN(I,J,R)-RMATRX(I,J,K) IF(DS PLAY. EQ. YSECTN) BINWIN (I , J, K)-RMATRX(I , K, J) IF(DSPLAY.EQ.XSECTN) BINWIN(I,J,K)=RMATRX(J,K,I) 603 CONTINUE IF(OUTLUN.NE. 60) WRITE (OUTLUN, 1020) DO 610 K-1,IRDIM RMlsK-l IF(FLOAT(RMI/6).NE.(FLOAT(KMl)/6.)) GOTO 604 IF(RMETHD. EQ. RMMRRT) WRITE (OUTLUN, 1070) IF(RMETHD.EQ.RMSMRT) WRITE(OUTLUN, 1080) IF(RMETHD. EQ. RMSSRT) WRITE (OUTLUN , 1090) 604 WRITE(OUTLUN,1OOO) DSPLAY, KMl IF(DSPLAY.EQ.ZSECTN) WRITE(OUTLUN,1O4O) + XSECTN, (L,L-O, 15) ,YSECTN IF(DSPLAY.EQ.YSECTN) WRITE(OUTLUN, 1040) + XSECTN,(L,L-0,15),ZSECTN IF(DSPLAY.EQ.XSECTN) WRITE(OUTLUN, 1040) + YSECTN,(L,L-O,15),ZSECTN DO 605 J-1,IRDIM JMl-J-l WRITE (OUTLUN, 101 O) JM1, (BINWIN (I, J,K) ,1-1 , IRDIM) 605 CONTINUE IF(FLOAT(KMI/3).NE.(FLOAT(KM1)/3.)) WRITE(OUTLUN,119O) 610 CONTINUE WRITE(OUTLUN, 1030) RTIME,DTIME,ET1ME WRITE (OUTLUN, 1020) GOTO 601 4:- * BINARY OUTPUT WINDOWING. g- 700 WRITE (60,1100) READ(6O, 1110) IANS IF(IANS.EQ.N) GOTO 900 71 WRITE(6O,1210) READ(6O,1110) DSPLAY WRITE(6O,1120) READ(6O,113O) THRESH DO 750 K-1,IRDIM DO 740 J=1,IRDDM DO 730 I-1,IRDIM BINWIN(I,J,R)-ATSIGN IF(DSPLAY.EQ.ZSECTN.AND.RMATRX(I,J,R).LT.THRESH) + BINWIN(I,J,R)sDOT IF(DSPLAY.EQ.XSECTN.AND.RMATRX(J,R,I).LT.THRESH) + BINWIN(I,J,R)=DOT IF(DSPLAY.EQ.YSECTN.AND.RMATRX(I,K,J).LT.THRESH) + BINWIN(I, J,R)=DOT 730 CONTINUE 740 CONTINUE 750 CONTINUE WRITE(6O, 1140) DSPLAY,LABEL(1),DSPLAY,LABEL(Z) ,DSPLAY,LABEL(3) , + (LABEL(I),I-1,16), + (LABEL(I),I-1,16),(LABEL(I),I-1,16) DO 800 J=I,I6 WRITE(60,1150) LABEL(J),(BINWIN(I,J,1),I-1,16), + LABEL(J),(BINWIN(I,J,2),I-1,16), + LABEL(J),(BINWIN(I,J,3),I-1,16) 800 CONTINUE WRITE(6O, 1180) WRITE(6O, 1140) DSPLAY,LABEL(7),DSPLAY,LABEL(S),DSPLAY,LABEL(9), + (LABEL(I),I-1,16),(LABEL(I),I-1,16),(LABEL(I),I-1,16) DO 810 J-1,16 WRITE(6O,1150) + LABEL(J),(BINWIN(I,J,7),I=1,16), + LABEL(J),(BINWIN(I,J,8),I-l,16), + LABEL(J),(BINWIN(I,J,9),I-1,16) 810 CONTINUE WRITE(6O,1180) WRITE(6O,114O) DSPLAY,LABEL(13),DSPLAY,LABEL(14),DSPLAY, LABEL(15) + ,(LABEL(I),I-1,16),(LABEL(I),I=1,16),(LABEL(I),I-1,16) DO 820 J-l,16 WRITE(6O,1150) LABEL(J),(BINWIN(I,J,13),I-1,16), + LABEL(J),(BINWIN(I,J,14),I-1,16), + LABEL(J),(BINWIN(I,J,15),I=l,16) 820 CONTINUE WRITE(60,1190) WRITE(60,1140) DSPLAY,LABEL(4),DSPLAY,LABEL(S),DSPLAY,LABEL(G), + (LABEL(I),I-l,l6), + (LABEL(I),I=1,16), (LABEL(I),I-l,l6) DO 830 J-1,16 WRITE(60,1150) + LABEL(J),(BINWIN(I,J,4),I-l,l6), + LABEL(J),(BINWIN(I,J,5),I-l,l6), + LABEL(J),(BINWIN(I,J,6),I-l,l6) 830 840 850 900 72 CONTINUE WRITE (60, 1180) WRITE (60, 1140) DSPLAY,LABEL(10) ,DSPLAY,LABEL(I l ) ,DSPLAY,LABEL(IZ + (LABEL(I),I-l , l6),(LABEL(I),I-1 , l6) , (LABEL(I ) ,I-l , 16) DO 840 J81,l6 WRITE(60,1150) LABEL(J),(BINWIN(I,J, lO),I-1,l6), + LABEL(J),(BINWIN(I,J, ll),I-1,16), + LABEL(J),(BINWIN(I,J,12),I-1,16) CONTINUE WRITE (60, 1180) WRITE(6O,116O) DSPLAY, LABEL(16),RMETHD,THRESH,(LABEL(I),I=1,16) DO 850 J-1,16 WRITE(6O,117O) LABEL(J),(BINWIN(I,J,16),I-1,16) CONTINUE WRITE(6O,1180) GOTO 700 CONTINUE END SUBROUTINE VOXGEN(DEPTH, VOXNUM,XM,YM,ZM) INTEGER VOXNUM C(MMON /RECON/ RMATRX(16,16,16), IRDIM, KSLAB GENERATES MS VOXEL COORDINATES ACCORDING TO CONVENTIONS OF SLAB GEOMETRY. MY- (VOXNUM-l ) /IRDIM MX- (VOXNUM-l )-MY*IRDIM RCENTR-FLOAT (IRDIM) /2 . -. 5 XM-MX-RC ENTR YM-MY-RC ENTR ZM'IDEPI'H RETURN END a-a-a-a-a-r- 8- 73 SUBROUTINE TRNSFM(VOXELX,VOXELY,VOXELZ,DAUG) REAL DAUG(3,4) CALLED WITH (VOXELX,VOXELY,VOXELZ) B (XM,YM,ZM); PERFORMS SPATIAL TRANSFORMATION SPECIFIED BY DAUG; RETURNS WITH (VOXELX,VOXELY,VOXELZ) - (X,Y,Z). fl- 1- X=VOXELX Y=VOXELY Z=VOXELZ VOXELXfiDAUG(1,l)*X +-DAUG(1,2)*Y +-DAUG(1,3)*Z +-DAUG(1,4) VOXELYaDAUG(2,1)*x +-DAUG(2,2)*Y +-DAUG(2,3)*Z +-DAUG(2,4) VOXELZ-DAUG(3,1)*X + DAUG(3,2)*Y +-DAUG(3,3)*Z +-DAUG(3,4) RETURN END SUBROUTINE MRRT(Ix,IY,Iz,VVALUE) COMMON /RECON/ RMATRX(16,16,16), IRDIM, KSLAB IMPLEMENTS MAXIMUM-REPLACEMENT RECONSTRUCTION TECHNIQUE. RMATRX(IX,IY,IZ)=AMAX1(RMATRX(Ix,IY,IZ),VVALUE) RETURN END SUBROUTINE DIRCOS(MPCODE,ANGLE,DAUG) REAL DAUG(3,4) INTEGER V,F,B,R,L DATA v,F,B,R,L/1HV,1HF,1HB,1HR,1HL/ DATA PI/3.141592654/ CALLED WITH TRANSDUCER.MEASUREMENT PLACEMENT CODE, ANGLE OF AIM WITH RESPECT TO VERTICAL, LOCATION IN RS (XO,Y0,ZO); RETURNS DIRECTION COSINE MATRIX.FOR MS-RS TRANSFORMATION. iii-”8'8- a. l INITIALIZE DAUG'TO COSINE OF PI/2. 100 PID2-PI/2. DO 100 I-1,3 DO 100 J-1,3 DAUG(I,J)-0. CONTINUE IF(MPCODE.NE.V) GOTO 200 74 CONSTRUCT MATRIX FOR "VERTICAL" MEASUREMENTS. DO 150 I=l,3 DAUG(I,I)-l. CONTINUE GOTO 900 IF(MPCODE.NE.F) GOTO 300 CONSTRUCT MATRIX FOR "FRONT" MEASUREMENTS. DAUG(1,1)-I. DAUG(2,2)-COS(ANGLE) DAUG(3,3)=DAUG(2,2) DAUG(2,3)-COS(PID2+ANGLE) DAUG (3 , 2)-=COS (PID 2-ANGLE) GOTO 900 IF(MPCODE.NE.B) GOTO 400 CONSTRUCT MATRIX FOR "BACK" MEASUREMENTS. 400 DAUG(1,1)--1. DAUG(2,2)-COS(PI:ANGLE) DAUG(3,3)-COS(ANGLE) DAUG(2,3)=COS(PID22ANGLE) DAUG(3,2)-DAUG(2,3) GOTO 900 IF(MPCODE.NE.R) GOTO 500 CONSTRUCT MATRIX FOR "RIGHT" MEASUREMENTS. DAUG (1 , 2)=-COS (ANGLE) DAUG (1 , 3)=COS (PID2+ANGLE) DAUG (2 , 1)--1 . DAUG (3 , 2)-COS (PID 2-ANGLE) DAUG (3 , 3)-DAUG (1 , 2) GOTO 900 IF(MPCODE.NE.L) STOP * CONSTRUCT MATRIX FOR "LEFT" MEASUREMENTS. DAUG (1 , 2)-COS (PI-ANGLE) DAUG (1 , 3)=COS (PID 2-ANGLE) DAUG (2, 1)-1. DAUG (3 , 2)-COS (PID 2-ANGLE) DAUG (3 , 3)-COS (ANGLE) 900 RETURN END 75 SUBROUT INE SMRT (IX, IY, IZ , VVALUE) COMMON IRECON/ RMATRX(16,16,16), IRDIM, KSLAB IMPLEMENTS SUMMATION RECONSTRUCTION TECHNIQUE. IF(RMATRX(IX, IY, IZ) .EQ. (-1.)) RMATRX(IX, IY, IZ)=0. RMATRX(IX, IY, IZ)-RMATRX(IX, IY, IZ)+VVALUE RETURN END SUBROUT INE SSRT (IX, IY, IZ , VVALUE) COMMON IRECON/ RMATRX(16,16,16), IRDIM, KSLAB IMPLEMENTS SUM OF SQUARES RECONSTRUCTION TECHNIQUE. IF(RMATRX(IX, IY, IZ) .EQ. (-1. )) RMATRX(IX, IY, IZ)'O. RMATRX(I X, IY, IZ)-RMATRX(I X, IY, IZ)+VVALUE **2 . RETURN END REFERENCES 10. 11. 12. 13. REFERENCES L. Filipczynski, A. Pawlowski and J. Wehr, Ultrasonic Methods of Testing Materials. Butterworths, London, 1966. W.J. McGonnagle, Nondestructive Testing. Second edition, Gordon and Breach, New York, 1961. J.C. Gore and S. Leeman, "Ultrasonic backscattering from human tissue, A realistic model". Phys. Med. Biol., Vol. 22, No. 2, 8.0. Aks and D.J. Vezzetti, "Improvements in the analysis of ultrasonic scattering from inhomogeneities in tissue". Ultrasonic Imaging, Vol. 1, No. 1, pp. 44-55, 1979. 8.0. Aks, remarks made during tutorial session at Fourth International Symposium on Ultrasonic Imaging and Tissue Characterization, National Bureau of Standards, Gaithersburg, MD, June 1979. J. Krautkramer, Ultrasonic Testing of Materials. Springer-Verlag, New York, 1969. E.E. Aldridge, "Ultrasonic holography". Engineering, Vol. 208, pp. 71-74, 19690 J.N. Butters, Holography and Its Technology. Peter Peregrinus, Ltd., London, 1971. J.L. Kreuger, "Ultrasonic holography for nondestructive testing". Materials Evaluation, Vol. 26, pp. 197-203, 1968. Y. Aoki, "Image reconstruction by computer in acoustical holo- graphy". Acoustical Holography, Vol. 5, pp. 551—572, 1974. S.J. Norton and M. Linzer, "Ultrasonic reflectivity tomography, Reconstruction with circular transducer arrays". Ultrasonic Imaging, Vol. 1, No. 2, pp. 154-184, 1979. S.J. Norton and M. Linzer, "Ultrasonic reflectivity imaging in three dimensions, Reconstruction with spherical transducer arrays". Ultrasonic Imaging, Vol. 1, No. 3, pp. 210—231, 1979. M. Onoe, M. Takagi, T. Masumoto and N. Hamano, "Graphic display for ultrasonic nondestructive testing". Acoustical Holography, V01. 4, pp. 299-315, 1972. 76 14. 15. 16. 17. 18. 19. 20. 77 J.B. Marion, Classical Dynamics of Particles and Systems. Academic Press, New York, 1971. Varian 74 System Handbook, Varian Data Machines, Inc., Irvine, CA, Feb. 1975. E.L. Madsen, J.A. Zagzebski and T. Ghilardi—Netto, "An anthropo— morphic torso section phantom for ultrasonic imaging". Med. Physics, Vol. 7, No. 1, pp. 43-50, 1980. D.E. Robinson, "Pulse-echo reconstruction techniques". Fifth International Symposium on Ultrasonic Imaging and Tissue Characterization, National Bureau of Standards, Gaithersburg, MD, June 1980. G.T. Herman and H.K. Liu, "Computer display of organ surfaces". A Review of Information Processing in Medical Imaging, ORNL/ BCTIC-Z, pp. 283-301, Biomedical Computing Technology Informa- tion Center, Oak Ridge National Laboratory, Oak Ridge, TN, 1977. G.T. Herman and H.K. Liu, "Display of three-dimensional informa- tion in computed tomography". J. Computer Assisted Tomography, V01. 1, No. 1, pp. 155-160, 1977. H.K. Liu, "Two- and three-dimensional boundary detection". Computer Graphics and Image Processing, Vol. 6, pp. 123-134, 1977.