Citation
Tomographic two-phase flow measurement using Compton scattering of gamma rays

Material Information

Title:
Tomographic two-phase flow measurement using Compton scattering of gamma rays
Creator:
Bodette, David E. ( Dissertant )
Jacobs, Alan M. ( Thesis advisor )
Carroll, Edward E. ( Reviewer )
Dugan, Edward T. ( Reviewer )
Hanrahan, Robert J. ( Reviewer )
Roan, Vernon P. ( Reviewer )
Place of Publication:
Gainesville, Fla.
Publisher:
University of Florida
Publication Date:
Copyright Date:
1986
Language:
English

Subjects

Subjects / Keywords:
Annular flow ( jstor )
Density ( jstor )
Gamma ray spectrum ( jstor )
Gamma rays ( jstor )
Inverted spectra ( jstor )
Noise spectra ( jstor )
Photons ( jstor )
Pixels ( jstor )
Spectral energy distribution ( jstor )
Spectral reconnaissance ( jstor )
Compton effect ( lcsh )
Dissertations, Academic -- UF -- Nuclear Engineering Sciences
Gamma rays ( lcsh )
Nuclear Engineering Sciences thesis Ph.D
Tomography ( lcsh )
Two phase flow ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Abstract:
Using Compton scattered gamma-rays to measure local void fraction was first suggested by Kondic and Hahn in 1970. The Compton scattered gamma-ray densitometers they suggested employ a single, narrow beam source and either a well collimated detector or an uncollimated detector. The collimated detector configuration only gives the local void fraction measurement in the small volume where the source and detector collimators intersect. The uncollimated detector configuration is a more efficient design since the local void fraction along the source beam’s path is measured in a single reading. A logical extension of the technique is to use wide beam illumination and uncollimated detectors to sample an even greater portion of the flow cross section. This report investigates and demonstrates several methods of inferring two-phase flow parameters using wide beam illumination coupled with two detectors placed symmetrically about he source and pipe. The spatial distribution of the fluid in a slice of the pipe is encoded with respect to energy in the singly scatter photon flux. Two basic techniques are detailed for decoding the spatial information: the method of spectral moments and the method of computed tomography (CT). Examination of the low-order moments of the singly scattered photon spectra from the two detectors provides sufficient information for flow regime identification. Flow asymmetries are revealed by comparison of the spectral moments of the two measured spectra. Further classification of the flow pattern is made on the basis of the first and second moments of the spectra. The real focus of this report is the adaptation and successful demonstration of CT techniques with Compton scattering. Three series expansion techniques are considered: the algebraic reconstruction technique (ART), the simultaneous iterative reconstruction technique (SIRT), and the iterative least squares (ILS) technique. Application of the modified ART, SIRT, and ILS algorithms to a hot-spot and a cold-spot reconstruction problem indicates that SIRT and ILS are more accurate than ART. A more extensive testing of the ILS algorithm for a variety of model flow regimes demonstrates the potential of Compton scatter tomography as a quantitative measurement technique.
General Note:
Typescript.
General Note:
Vita.
Thesis:
Thesis (Ph.D.)--University of Florida, 1986.
Bibliography:
Bibliographies : leaves 191-195.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
030321815 ( alephbibnum )
16656185 ( oclc )

Downloads

This item has the following downloads:

Binder1 ( .pdf )

UF00090205_00001.pdf

00006.txt

00199.txt

00206.txt

00026.txt

00047.txt

00080.txt

00058.txt

00105.txt

00060.txt

00054.txt

00092.txt

00051.txt

00177.txt

00055.txt

00061.txt

00153.txt

00162.txt

00137.txt

00205.txt

00183.txt

00067.txt

00142.txt

00181.txt

00037.txt

00033.txt

00100.txt

00096.txt

00145.txt

Binder1_pdf.txt

00108.txt

00174.txt

00062.txt

00002.txt

00112.txt

00146.txt

00076.txt

00057.txt

00148.txt

00182.txt

00158.txt

00087.txt

00066.txt

00186.txt

00073.txt

00075.txt

00194.txt

00007.txt

00127.txt

00027.txt

00063.txt

00114.txt

00091.txt

00071.txt

00120.txt

00059.txt

00136.txt

00150.txt

00042.txt

00012.txt

00201.txt

00156.txt

00125.txt

00023.txt

00167.txt

00039.txt

00122.txt

00163.txt

00133.txt

00210.txt

00072.txt

00081.txt

00020.txt

00038.txt

00188.txt

00179.txt

00193.txt

00151.txt

00101.txt

00011.txt

00190.txt

00160.txt

00034.txt

00010.txt

00083.txt

00157.txt

00143.txt

00024.txt

00110.txt

00093.txt

00117.txt

00152.txt

00184.txt

00022.txt

00204.txt

00119.txt

00189.txt

00168.txt

00111.txt

00154.txt

00207.txt

00019.txt

00203.txt

00126.txt

00135.txt

00172.txt

00191.txt

00170.txt

00169.txt

00070.txt

00032.txt

00138.txt

00068.txt

00107.txt

00128.txt

00140.txt

00064.txt

00008.txt

00035.txt

00095.txt

00200.txt

00090.txt

00196.txt

00016.txt

00116.txt

00118.txt

00005.txt

00103.txt

00208.txt

00166.txt

00197.txt

00017.txt

00139.txt

00178.txt

00097.txt

00050.txt

00121.txt

00085.txt

00195.txt

00018.txt

00098.txt

Copyright.txt

00209.txt

00113.txt

00052.txt

00144.txt

00084.txt

00069.txt

00134.txt

00004.txt

00088.txt

00187.txt

00029.txt

00175.txt

00074.txt

00132.txt

00077.txt

00041.txt

00053.txt

00164.txt

00198.txt

00104.txt

00185.txt

00115.txt

00078.txt

00149.txt

00141.txt

00131.txt

00021.txt

00028.txt

00031.txt

00009.txt

00046.txt

00147.txt

00044.txt

00013.txt

00001.txt

00109.txt

00099.txt

00102.txt

00180.txt

00040.txt

00129.txt

00094.txt

00159.txt

00014.txt

00086.txt

00130.txt

00049.txt

00079.txt

00048.txt

00165.txt

00123.txt

00065.txt

00106.txt

00015.txt

00056.txt

00192.txt

00045.txt

00161.txt

00171.txt

00176.txt

00173.txt

00202.txt

00030.txt

UF00090205_00001_pdf.txt

00089.txt

00082.txt

00155.txt

00036.txt

00124.txt

00043.txt

00025.txt

00003.txt


Full Text














TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT
USING COMPTON SCATTERING OF GAMMA RAYS






By

DAVID E. BODETTE


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY


UNIVERSITY OF FLORIDA


1986



































Copyright 1986

by

David E. Bodette









ACKNOWLEDGMENTS

I thank Dr. Alan M. Jacobs, my adviser, for his patient

guidance and careful attention to detail.

This work was partially supported by a Liquid Metal

Fast Breeder Reactor Technology Research Award from the

Argonne National Laboratory under contract with the United

States Department of Energy (DOE). The author was also

supported by the DOE Energy Graduate Traineeship Program.


iii










TABLE OF CONTENTS


ACKNOWLEDGMENTS . . . . . . . . . .

LIST OF TABLES . . . . . . . . . .

LIST OF FIGURES . . . . . . . . .

ABSTRACT . . . . . . . . . . . .

CHAPTERS

1 INTRODUCTION . . . . . . . .

2 THE FORWARD PROBLEM . . . . . . .

Collimated Point Source/Collimated Detector .
Collimated Point Source/Uncollimated Detector
Uncollimated Point Source/Uncollimated Detector
Parallel Beam Source/Uncollimated Detector. .


Discrete Forward Equations


3 METHOD OF SPECTRAL MOMENTS. . . .

Description of Method . . . . .
Procedure and Results . . . . .
Interpreting the Moments. . . .
A Detailed Example . . . .
Using the Moments . . . . .

4 THE RECONSTRUCTION PROBLEM. . . .

General Solution Methods . . . .
A Closer Look . . . . . . .
Backprojection . . . . . .
Monte-Carlo Backprojection . . .
Matrix Methods . . . .
Algebraic Reconstruction Technique
Simultaneous Iterative
Reconstruction Technique
Iterative Least-Squares Technique .
Comparison of Reconstruction Techniques

4 COMPTON SCATTER TOMOGRAPHY . . .

An Example . . . . . . .
Accuracy . . . . . . .
Resolution . . . . . . .
Image Contrast . . . . . .


. . 35

. . 35
. . 42
. . 44
. . 46
. . 52

. . 54

. . 56
. . 58
. . 59
. . 64
. . 69
. . 72

. . 76
. . 77
. . 79

. . 90

. . 90
. . 92
. . 106
. . 109


Page

iii

vi

vii

xi




1

10

11
15
17
28


. . . . . 29











6 SUMMARY AND CONCLUSIONS . . . . .

Summary . . . . . . . . .
A Practical Application . . . . .
Conclusions and Future Research . . .

APPENDICES

I MONTE-CARLO TRANSPORT PROGRAM . . . .

II INSIGHTS AND FINER POINTS . . . . .

Parallel Versus Diverging Source Beams. .
Multiple Scatter Versus Detector Position .
Source Placement . . . . . .
Convergence Criterion . . . . . .


III RECONSTRUCTION COMPUTER PROGRAMS ..

Running the Programs . . . . . .
Flow of Programs . . . . . . .

IV MAXIMUM ENTROPY METHOD . . . . .

V DATA . . . . . . . . . .

REFERENCES . . . . . . . . . .

BIBLIOGRAPHY . . . . . . . . . .

BIOGRAPHICAL SKETCH . . . . . . . .


Page

118

119
121
126



130

133

133
137
142
143


S. 144

S. 144
. 145

S. 146

S. 149

S. 191

S. 195

S. 196









LIST OF TABLES


Page

2.1 Multiply scattered photon spectra . . .. 22

3.1 Spectral moment results . . .. . . 45

4.1 Data for hot-spot and cold-spot tests ... . 83

5.1 Test parameters and values . . .. . . 98

5.2 Summary of reconstruction results. . . . 105

5.3 Parameters affecting image resolution . . 107









LIST OF FIGURES


Page
2.1 Collimated source/collimated detector
configuration 12

2.2 Collimated source/uncollimated detector
configuration 16

2.3 Uncollimated source/uncollimated detector
configuration 19

2.4 Illustration of scattering angle blur 25

2.5 Illustration of an isogonic segment 27

2.6 Parallel beam source/uncollimated detector
configuration 30

2.7 Overlaying of pipe cross section with an n x n
pixel grid 30

2.8 Feathering of isogonic segments 34

3.1 Compton scatter densitometry configuration 36

3.2 Spectral moments of a single-scattered spectrum 39

3.3 Source-pipe-detector configuration 44

3.4 Illumination of stratified flow regime 44

3.5 Pipe illumination for various beam widths 48

3.6 Spectrum of upper detector versus beam width 49

3.7 Spectrum of lower detector versus beam width 50

4.1 Conventional CT scanner 60

4.2 Monte-Carlo backprojection 67

4.3 Spectra from Monte-Carlo backprojection 68

4.4 Flow chart of reconstruction algorithms 71

4.5 Example of an ART solution for two pixel image 75

4.6 Hot-spot and cold-spot test configuration 80

4.7 Hot-spot spectra 81

4.8 Cold-spot spectra 82


vii









Page
4.9 Hot-spot reconstruction 84

4.10 Cold-spot reconstruction 85

4.11 Convergence of reconstruction algorithms 87

5.1 Calculated spectrum using incorrect source
strength value 91

5.2 Reconstruction using incorrect source
strength value 93

5.3 Experimental setup 95

5.4 Flow regime models 99

5.5 Singly scattered photon spectrum 100

5.6 Reconstruction with and without noise 101

5.7 Spectra with noise added 103

5.8 Detected spectrum with pipe wall effects 117

6.1 Spectrum from a 4 inch pipe 122

6.2 Detail of spectrum from the pipe 123

6.3 Effect of signal-to-noise ratio 125

6.4 Detector response versus bubble size 127

I.1 Example geometry of photon trajectory 132

II.1 Detector response for diverging source beam 135

11.2 Detector response for parallel source beam 136

11.3 Detector configuration for multiple scatter
versus detector position 138

11.4 Spectra for multiple scatter versus detector
position 139

V.1 Spectra from annular flow 151

V.2 Spectra from annular flow 152

V.3 Spectra from annular flow 153

V.4 Spectra from inverted annular flow 154


viii









Page
V.5 Spectra from inverted annular flow 155

V.6 Spectra from inverted annular flow 156

V.7 Spectra from bubbly flow 157

V.8 Spectra from bubbly flow 158

V.9 Spectra from bubbly flow 159

V.10 Spectra from stratified symmetric flow 160

V.11 Spectra from stratified symmetric flow 161

V.12 Spectra from stratified symmetric flow 162

V.13 Spectra from stratified asymmetric flow 163

V.14 Spectra from stratified asymmetric flow 164

V.15 Spectra from stratified asymmetric flow 165

V.16 Annular flow pattern 166

V.17 Reconstruction of annular flow 167

V.18 Reconstruction of annular flow 168

V.19 Reconstruction of annular flow 169

V.20 Inverted annular flow pattern 170

V.21 Reconstruction of inverted annular flow 171

V.22 Bubbly flow pattern 172

V.23 Reconstruction of bubbly flow 173

V.24 Reconstruction of bubbly flow 174

V.25 Reconstruction of bubbly flow 175

V.26 Stratified symmetric flow pattern 176

V.27 Reconstruction of stratified symmetric flow 177

V.28 Stratified asymmetric flow pattern 178

V.29 Reconstruction stratified asymmetric flow 179

V.30 Example of smooth convergence 180









Page
V.31 Example of divergence 181

V.32 Example of divergence 182

V.33 Example of oscillating divergence 183

V.34 Example of non-convergence 184

V.35 Setup for laboratory experiments 185

V.36 Experimental spectrum of annular flow 186

V.37 Experimental and Monte-Carlo spectra 188

V.38 Experimental spectrum of stratified flow 189

V.39 Experimental spectrum of bubbly flow 190















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy


TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT
USING COMPTON SCATTERING OF GAMMA RAYS

By

DAVID E. BODETTE

December, 1986

Chairman: Dr. Alan M. Jacobs
Major Department: Nuclear Engineering Sciences

Using Compton scattered gamma-rays to measure local

void fraction was first suggested by Kondic and Hahn in

1970. The Compton scattered gamma-ray densitometers they

suggested employ a single, narrow beam source and either a

well collimated detector or an uncollimated detector. The

collimated detector configuration only gives the local void

fraction measurement in the small volume where the source

and detector collimators intersect. The uncollimated

detector configuration is a more efficient design since the

local void fraction along the source beam's path is measured

in a single reading. A logical extension of the technique

is to use wide beam illumination and uncollimated detectors

to sample an even greater portion of the flow cross section.

This report investigates and demonstrates several

methods of inferring two-phase flow parameters using wide

beam illumination coupled with two detectors placed










symmetrically about the source and pipe. The spatial

distribution of the fluid in a slice of the pipe is encoded

with respect to energy in the singly scattered photon flux.

Two basic techniques are detailed for decoding the spatial

information: the method of spectral moments and the method

of computed tomography (CT).

Examination of the low-order moments of the singly

scattered photon spectra from the two detectors provides

sufficient information for flow regime identification. Flow

asymmetries are revealed by comparison of the spectral

moments of the two measured spectra. Further classification

of the flow pattern is made on the basis of the first and

second moments of the spectra.

The real focus of this report is the adaptation and

successful demonstration of CT techniques with Compton

scattering. Three series expansion techniques are

considered: the algebraic reconstruction technique (ART),

the simultaneous iterative reconstruction technique (SIRT),

and the iterative least squares (ILS) technique.

Application of the modified ART, SIRT, and ILS algorithms to

a hot-spot and a cold-spot reconstruction problem indicates

that SIRT and ILS are more accurate than ART. A more exten-

sive testing of the ILS algorithm for a variety of model

flow regimes demonstrates the potential of Compton scatter

tomography as a quantitative measurement technique.


xii














CHAPTER 1
INTRODUCTION



The case for developing improved non-intrusive measur-

ing technologies for sensing material/phase distribution in

two-phase flow has been well argued by many authors. The

problem is one which spans a broad number of fields inclu-

ding manufacturing, engineering, security, and medicine.

Each potential application imposes its own peculiar demands,

which are often at odds with constraints in other applica-

tions. This situation has resulted in many diverse

techniques being pursued and often optimized for each

specific problem. Examples of the technologies employed

include ultrasonics, microwaves, lasers, optics, and irra-

diation with particles, X rays, or gamma rays.

One specific problem encountered in the nuclear power

industry and its associated research industry is to measure

the void fraction and/or liquid distribution (flow regime)

of two-phase flows inside steel walled pipes. Any intrusive

instrument which disturbs the flow pattern is certainly not

desirable, but an even more strict requirement often encoun-

tered in these circumstances is that penetrations through

the pipe wall are not permitted. Another requirement of any

measuring instrument is that its output should be unambi-

guous. In the case of two-phase flow void fraction









measurement this means the entire flow cross section must be

interrogated. There may also be physical limitations in

terms of size and orientation, as well as access to the

pipe.

A perusal of the recent literature on measurements in

two-phase flow [1, 2, 3, 4] reveals that the choice of

techniques for making non-intrusive local void fraction

and/or phase distribution measurements is limited to irra-

diation with particles, gamma rays, and laser light. When

intrusion is not an issue, a variety of probe techniques are

available for local property measurements. Reference [4]

lists four basic probe methods which appear to be reliable

enough for general use: electrical probes, optical probes,

thermal anemometers, and micro-thermocouples. Laser

scattering and interferometry is suggested as a very

promising method when the flow boundary is not opaque.

There are techniques for making average void fraction

measurements which might be adaptable to the local void

fraction problem. Ultrasonics is capable of measuring the

average void fraction in very thick walled pipes, but only

for low void fractions (i.e. < 0.20) [3]. An array of

ultrasonic sensors around a pipe might produce an image of

the flow provided the voids are not too big or too many.

Microwave resonant cavities are also capable of measuring

average void fraction. Using microwaves to sense local void

fraction would probably run into difficulties from the metal

pipes [1]. Nuclear magnetic resonance (NMR) is currently










receiving a lot of attention in the medical field and

certainly has potential for two-phase flow. The medical NMR

machines involve a radio-frequency field and a large

magnetic field which cause problems in nearby electronic

equipment. The metal pipes seem to be another source of

potential problems in implementing NMR for local

measurements; however, plastic pipes can be used to avoid

such problems [5].

Gamma-ray attenuation is the predominant method of

choice for making average void fraction measurements and is

often used as the calibration standard for other methods.

Gamma-ray attenuation instrumentation and techniques are

recognized as having reached the highest level of

development. The reason for this dominance is that the

method places a minimum of constraints on the user. The

method is non-intrusive, fairly insensitive to flow pattern,

and functions over a broad range of void fractions, flow

regimes, and temperatures.

There are many ways of utilizing gamma rays or neutrons

to measure local void fraction and phase distribution. The

basic methods of implementing tomographic two-phase flow

measurement with gamma rays or neutrons can be grouped into

the three following categories:


1) Measure the emission of gamma rays or neutrons from-
activation products in the fluid, induced
excitation, or injected radiotracers.

2) Measure the attenuation of narrow or wide source
beams in traversing the pipe.










3) Measure the singly scattered gamma rays or neutrons
from a primary source beam.

Methods from category 1 could be roughly classified as

emission computed tomography (ECT) techniques. An example

of a category 1 method is a single photon ECT system which

examines the ~6 MeV gamma rays from the activation-decay

scheme pair



16(n,p)N16
and
N16 O16 + y + 8.



Although recent advances in the analytical techniques of ECT

and in instrumentation suggest this as a viable method, it

appears that only cursory attention has been given to the

technique [6]. An ECT system could also be employed using

injected radiotracers; however, the quantity of radio-

nuclides required poses a serious drawback [7].

The most familiar tomographic instruments comprise

category 2 and are generally referred to as computed axial

tomography (CAT) scanners. The basic hardware setup of CAT

scanners involves placing a source on one side of the pipe

and a detector on the opposite side of the pipe. The whole

source-detector assembly rotates in a controlled manner

around the pipe while the attenuation of the source beam at

different orientations is recorded. The rate at which CAT

scanners can rotate about the pipe is a limiting factor.

CAT scanners of the newest generation do not employ any









moving parts, but instead use a magnetically focused flying

spot electron beam and a fixed circular detector array [8].

These fifth generation CAT scanners can acquire an image in

approximately 1/10 second and are very promising where cine

tomography is necessary. The capital investment in these

new CAT scanners is significant and combined with the fact

that they are not compact or portable means they are not a

viable alternative in many applications.

Category 3 techniques are closely related to ECT

except that the emitted radiation does not come from

radioactive material within the fluid, but from radiation

scattered out of an external illuminating source beam.

Gamma-ray and neutron sources are possible candidates, and

methods based on both are actively being developed [9].

The earliest reported efforts with scattered radiation

focused on gamma rays scattering and involved measuring the

local density in a small volume [10]. A well collimated

source illuminated an object and a well collimated detector

measured the Compton scattered gamma rays coming from the

small volume where the source and detector collimators'

fields-of-view intersected. The output of the detector was

proportional to the electron density in the small volume

which in turn is proportional to the material density for a

constant material composition. The Compton scattering

method of local density measurement developed, at first,

independently in the medical and industrial imaging fields.

Lale published two articles, one in 1959 [11] and another in








1968 [12], which demonstrated how internal organs could be

imaged by exploiting the Compton scattering of gamma rays.

In 1970, Kondic and Hahn [13] published the first suggested

engineering application of Compton scattering, that of

making local density measurements in two-phase flow. They

proposed both a collimated source/collimated detector (as

Lale) and a collimated source/uncollimated detector arrange-

ment. The uncollimated detector offered the improvement

that local density values were given for points along the

illuminated path of the source beam from a single

measurement. Acting independently, Farmer and Collins [14,

15] proposed and demonstrated the collimated source/uncolli-

mated detector arrangement in determining anatomical cross

sections. Their images were only qualitative since the

calibration procedures used failed to completely remove the

effects of multiple scattering or variations in attenuation.

The analytic method for making quantitative measurements

using this source-detector arrangement by deconvolving the

substantial multiple scatter blur was developed by Anghaie

[16, 17].

The use of Compton scattering to make local density

measurements in two-phase flows was pursued earnestly in the

late 70's and in to the 80's. Almost exclusively the focus

was on the collimated source and either collimated or uncol-

limated detector arrangement. Tomographic images produced

with either of these arrangements required scanning the

source and/or detector to cover all the points in the plane









of interest. Kondic [18] in 1978 first proposed using an

uncollimated or wide beam source coupled with an uncolli-

mated detector to image an entire flow cross section. Such

an arrangement made more efficient use of the source photons

than the previous source-detector arrangements and required

no moving parts. Both of these attributes made the arrange-

ment more adapted to the high framing rates necessary for

transient analysis. However, no critical study or viable

implementation of the wide source beam/uncollimated detector

arrangement has appeared in the literature.

Imaging with scattered neutrons was developed in the

late 70's and early 80's independently of and somewhat concur-

rently of the present work with scattered gamma rays. The

use of neutrons to make local void fraction measurements was

suggested by Banerjee and Lahey [19] and separately by

other researchers. Neutron scattering tomography was

developed by Hussein [20] as his Ph.D. thesis. He used a

wide parallel beam of 14 MeV neutrons to illuminate a

simulated two-phase flow and then measured the fluence and

energy of the scattered neutrons using small organic scinti-

lator detectors at different positions around the flow

volume.

The focus of the present dissertation is the develop-

ment and demonstration of the algorithms for reconstructing

the local density distribution of a two-phase flow using

Compton scattered gamma rays. The specific experimental

setup considered uses a wide-beam collimated gamma-ray









source to illuminate the flow and energy sensitive detectors

to examine the emerging singly Compton scattered photon

flux. To date, the reconstruction techniques proposed for

this setup have been based on deterministic algebraic solu-

tions and have not addressed the problems of multiple scat-

tering, statistical error, image blur, and other errors.

The contribution of this author is to apply the powerful

reconstruction algorithms developed in the medical and

optical imaging fields to this problem and to critically

assay the effects of__ multiple scattering and statistical

errors on the reconstruction.

Tomographic imaging with Compton scattered gamma rays

has two distinct features which distinguish it from

previously encountered tomographic reconstruction problems.

First, the imaging processes dealt with in the medical and

optic fields have been linear or could be linearized by an

analytic technique. Imaging the phase distribution in a

two-phase flow with Compton scattered gamma rays is a non-

linear problem. Second, the combination of using a point

source and the Compton scattering process introduces a

unique curved geometry to the flow reconstruction process.

The forward problem of predicting the emerging singly

Compton scattered photon flux given the local density dis-

tribution is discussed in Chapter 2. Chapter 3 details a

method of inferring tomographic information about a two-

phase flow by examining the spectral moments of the Compton

scattered photons [21]. Chapter 4 develops several more








sophisticated tomographic reconstruction techniques and

compares their performance under ideal conditions using data

generated by Monte-Carlo calculations. Chapter 5 evaluates

Compton scatter tomography under more realistic conditions

and discusses some of the considerations in designing a

system. Chapter 6 summarizes the significant results of this

study and offers directions for future work.















CHAPTER 2
THE FORWARD PROBLEM



Fundamental to the problem of reconstructing the phase

or density distribution using Compton scattered gamma rays

is being able to model the gamma ray interaction with the

density field. The aim of this chapter is to describe the

basic equations predicting the emerging scattered photon

flux from a known density distribution illuminated by a

point gamma-ray source. The basic concepts of using gamma-

ray scattering to make local density measurements are

introduced first with the collimated point source/collimated

detector arrangement. Once the equations for analyzing this

setup are developed, the analysis is extended to collimated

point source/uncollimated detector arrangements by introduc-

ing the equivalent approach of energy collimation. The

analysis is extended further to the wide-beam collimated

point source/uncollimated detector arrangement. Lastly,

parallel source beam illumination is considered and the

differences between it and the diverging wide-beam

illumination are discussed. Throughout, a detailed analysis

is only given to the singly scattered photon flux. The

presence of the multiply scattered photon flux is included

in all the equations and discussed in general terms.

However, for the case of the wide-beam collimated source









/uncollimated detector configuration a specific method is

given for removing the multiple scatter component.


Collimated Point Source/Collimated Detector

Figure 2.1 illustrates the arrangement of a collimated

source and collimated detector for making local density

measurements. Monoenergetic gamma rays leave the source and

are collimated to a thin pencil beam which illuminates a

chord of the pipe. A similar collimator on the detector is

aligned so that its field-of-view intersects the source

beam at the small volume denoted by V. The detector

senses the source gamma rays scattered in the volume V.

The detector response, D(E), is


D(E) = ESOPwAiAoAZ + M(E)


(2.1)


where S

p

pw


Az

Ai and Ao






M(E)


is the source activity;

is the density of the fluid in the volume V;

is the incoherent differential scattering cross
section;

is the path length of source photons in V;

are the total attenuation functions, including
interaction attenuation and the solid angles of
view, along the incoming and outgoing paths
respectively;

is the detector efficiency; and

represents the contribution to the spectrum of
multiply scattered photons.













Pipe


Collimator


Source


Detector


Figure 2.1: A collimated source/collimated detector gamma-
ray scattering configuration for measuring the
density in the volume V.









The functional form of both attenuation functions is



Ax = Cx*exp(-f pxp(z)dz) (2.1a)
zx



where Cx contains the geometric or field-of-view factors,

P, is the mass attenuation coefficient,

p(z) is the local density at z, and

x is a dummy subscript representing i and o for the
incoming and outgoing paths, respectively.

The line integral in equation (2.la) is carried out either

along the incoming path, ki, from the source to the volume V

or along the outgoing path, o from V to the detector. The

value of the integral is in units of optical path length.

The predominant interaction of gamma rays with energies

between 0.1 and 1.0 MeV in water is Compton scattering. The

energy of the singly Compton scattered photons arriving at

the detector is related to the polar scattering angle, e, by

the kinematic relation



E = Ei/(l + (Ei/mec2)(1 cos 8)) (2.2)



where Eo is the energy of the scattered photon,

Ei is the energy of the source photon, and

mec2 is the rest mass energy of an electron.


The singly scattered photons, represented by the first

expression on the right hand side (RHS) of equation (2.1),

all arrive at the detector after scattering through










approximately the same angle and thus have approximately the

same energy. Assuming the attenuation along the incoming

and outgoing paths remains constant in time, the detector

response to the singly scattered photons is proportional to

the density in the small volume V.

The multiply scattered photons arriving at the

detector, unlike the singly scattered photons, have a broad

range of energies since they can suffer any number of

collisions (with concomitant changes in energy) prior to

reaching the detector. The functional form of the multiply

scattered photon component of the detector response is



M(E)= ffff ('(r ,E', ')K(r',E',Q-;rd,E,Q)dr'dE'dQ'drdd ( (2.1b)



where M(..)is the angular photon flux and K(..) is the

transport kernel. The kernel can be viewed as the probabi-

lity that a photon at r going in direction n' with an

energy E' suffers a scattering interaction at r'; the scat-

tered photon goes in direction Q toward the detector with an

energy E; the scattered photon reaches point rd in the

detector's active volume, interacts with the detector and is

counted. The integral with respect to r is carried out

only over the detector collimator's field-of-view. The

detector collimator serves not only to pick out the volume

of interest, V, but to reduce the magnitude of the multiply

scattered photon spectral component by restricting the









region of space where the photons can suffer their last

collision prior to reaching the detector.


Collimated Point Source/Uncollimated Detector

The function of the detector collimator in picking out

the small volume V can also be achieved by electronic colli-

mation. An example is to use a high energy resolution

detector coupled to a multi-channel analyzer (MCA) instead

of the detector collimator. The point of scattering of a

singly scattered photon is related to its energy by equa-

tion (2.2) and the polar scattering angle as determined from

the measurement geometry. The singly scattered photons

arriving at the detector now have a range of energies. The

detector response becomes



D(E) = Sop(x)pwAi(x)Ao(x)e(E) + M(E) (2.3)



where S, p, 1,, Ai, Ao, and M are as defined for equa-

tion (2.1) and x is the point of scattering. The x-axis is

colinear with the source beam illumination path (Figure

2.2). Note that x and e can be expressed in terms of E

using equation (2.2) as


2
Yd
x = xd -
1 cos2e
(2.4)
cos e = 1 + (mec2/Ei)(l Ei/E)













Pipe


Collimator




Source- -


x- d


Detector


Figure 2.2: Configuration of a collimated source/uncollimated
detector gamma-ray scattering densitometer.








where xd and yd are the detector coordinates illustrated in

Figure 2.2. The first term on the RHS of equation (2.3) is

strictly a function of a single variable for a fixed

geometry.

The multiple scattering term in equation (2.3) has the

same functional form as equation (2.1) except that the inte-

gration with respect to r is over all space since there is

no detector collimator. Therefore, the magnitude of the

multiply scattered photon component of the spectrum is

greater with an uncollimated detector than with a collimated

detector. Offsetting this increased measurement structured

noise is the fact that the singly scattered spectrum now

contains the local density information for the entire source

beam illumination path. The information which would have

required many measurements with a collimated detector is

available from a single measurement with an uncollimated

detector.


Uncollimated Point Source/Uncollimated Detector

Just as eliminating the detector collimator increases

the information contained in a single measurement, so does

eliminating the narrow source collimator. An uncollimated

or a wide angle-of-view collimated source can illuminate the

entire flow cross section. The detector response then

contains the local density information from every point in

the flow cross section. The manner in which the local

density information is encoded in the detector response,









however, complicates the data processing necessary to

reconstruct the local density distribution.

Figure 2.3 illustrates a wide beam source Compton scat-

tering densitometer setup. The curves drawn across the pipe

represent lines of constant scattering angle, called

isogonic lines by Kondic [18]. Primary source photons

scattering at points along an isogonic line all scatter

through the same angle in order to reach the detector.

Using elementary geometry it is easy to show that the

isogonic lines are all segments of circles which pass

through both the source and the detector. The axis labeled

Y in Figure 2.3 is the loci of all the centers of these

circles. The radius, r, and the location of the center-

point, yc, of each isogonic circle as a function of the

scattering angle, 8, are



r = SD sin 8
2

(2.5)
y SD cot e
C- 2



where SD is the distance between the source and the

detector.

The detector response is now represented by an integral

equation



D(E)= f D(x,y)dk + M(E) (2.6)
Z(E)









Pipe


Source A


TOP VIEW


ELEVATION VIEW








Figure 2.3: Configuration of an uncollimated source/
uncollimated detector gamma-ray scattering
system for tomography.









where M(E) is the multiply scattered photon spectrum,

D(x,y) is the contribution to the detected photon
spectrum of singly Compton scattered photons
scattered at point (x,y), and

dZ is the differential line element along the
isogonic contour (E).

The first term on the RHS of equation (2.6) is essentially

the integral of the first term on the RHS of equation (2.3).

The second term, as before, represents the multiply

scattered photon contribution to the spectrum and has the

same functional form as equation (2.1b).

The multiple scatter component can comprise approxi-

mately one-third of the total signal [22, 23] and must be

removed. It may be possible to remove the multiple scatter-

ing component analytically since the transport kernel is a

known function and the angular photon flux can be generated

from an order-of-scattering calculation once a good estimate

of the density distribution is obtained. However, this

would be a cumbersome process and does not appear to be

necessary as there appears to be an acceptable alternative

approximation.

A multiply scattered photon arriving at the detector

with energy E could have come from any point within a large

volume of the pipe. The structure of the fluid distribution

imposed on the multiply scattered photon spectrum would tend

to be averaged out because of the many possible paths that

the photon could have taken prior to arriving at the

detector with energy E. The result is that for wide beam

source/uncollimated detector arrangements the multiply









scattered photon spectrum should be smooth and a simple

background subtraction procedure can be employed to remove

it.

The validity of this intuitive argument has been

tested by a series of Monte-Carlo calculations (see Appendix

I). The multiply scattered photon spectrums have been com-

puted for five flow regimes at both a high and low void

fraction. The calculations assume that the flow is

contained in a pipe whose thickness and composition are such

that it can be neglected. The results, summarized

pictorially in Table 2.1, show that to a first order approx-

imation equation (2.6) can be rewritten as



D(E)= J D(x,y)dk + B(E). (2.6*)
R(E)



where B(E) is some function of energy which is fitted to
the smooth multiply scattered photon
spectrum, and

D(x,y) is the contribution to the detected photon
spectrum of singly Compton scattered photons
scattered at point (x,y).

The function D(x,y) can be written more explicitly as



D(x,y)= SO (E)pw(Ei,Eo)p(x,y)Ai(x,y)Ao(x,y) (2.7)



where e(E) is the detector efficiency,

P, is the incoherent differential scattering cross
section for water,

p(x,y) is the local fluid density at (x,y),








TABLE 2.1

Monte-Carlo multiple scatter calculation results.


Stratified


Bubbly


Annular


Inverted Annular









Table 2.1 continued.









- - multiple
scatter

single
scatter ,"








Ai and Ao are as defined for equation (2.1), and

S is the source strength.


There are several implicit assumptions expressed in

equations (2.6 ) and (2.7). They are that a) the

divergence and width of the source beam in the direction

perpendicular to the imaging plane is small, b) the source

is a point body, c) the detector is a point body, and d) the

effects of the detector electronics and MCA are neglected.

Careful construction of the source collimator can insure

that a) is true. Including the effects of b), c), and d),

poses no serious analytic hurdles as far as the forward

problem of predicting the detector output given the density

distribution. However, these details have serious

implications on the resolution obtainable when trying to

reconstruct the density distribution from the detector

output.

,Equation (2.4) provides a one-to-one correspondence

between the energy of a detected singly scattered photon and

the point of scattering. The relationship is true only if

the source and detector are point bodies. If the so-uce ana

detector have finite extent, then there is some uncertainty

associated with the scattering angle as depicted in

Figure 2.4. This means that photons scattering at a point

in the density field do not arrive at the detector with a

single energy but with an energy over a small range,


AE = (@E/@e)Ae















so






e d
































Figure 2.4: Illustration of the scattering angle blur
problem arising from the use of non-point
source and detector.









This "blurring" effect is analogous to taking a picture with

an out-of-focus camera. The amount of uncertainty in the

scattering angle is not constant over the entire density

field nor even along an isogonic line. Consequently, the

blurring cannot be expressed as the convolution of some

spread function with the spectrum and cannot De removed by

a straignt forward deconvolutional procedure.

There is also blurring of the detected photon spectrum

by the detector and electronics. Unlike the blurring

discussed above, this blurring is generally energy invariant

and can be removed by a deconvolutional procedure. Methods

for removing the detector system blur and correcting for the

energy dependent efficiency are well developed and it is

assumed that the measured data are automatically corrected

for these effects.

Another aspect of the data collection process which

affects image quality is the use of a multi-channel analyzer

(MCA) to record the photon spectrum. The impact of collec-

ting the spectrum data in an MCA is that the spectrum is not

available as a continuous function but as a discrete set of

values. That is, D(E) becomes the set {Dc;c=l...I} where



Dc = D(E)dE.
AEc



It is now convenient to introduce the isogonic segment as

shown in Figure 2.5. An isogonic segment is the region

bounded by the isogonic curves which correspond to the

























































Figure 2.5: Illustration of an isogonic segment. (Reproduced
by permission. Copyright ISA 1985 [21].)










energy limits of a single channel in the MCA. Any features

of the spectrum smaller than the width of a single MCA

channel are smoothed over and cannot be recovered.

Generally this does not pose a problem since modern equip-

ment can collect data with very narrow channel widths.


Parallel Beam Source/Uncollimated Detector

The analysis so far has dealt exclusively with a small

source and a diverging beam. It is relevant to extend the

analysis to parallel source beams since work has been done

using parallel beam illumination with neutron scattering

[20]. Figure 2.6 illustrates one possible way of

constructing a parallel beam source from an array of

narrow beam sources. Alternatively, a parallel beam

source can be seen as the result of moving a wide beam

point source infinitely far away from the pipe. As the

source is moved away, the source-to-detector distance goes

to infinity and by equation (2.5) the radius of an

isogonic curve also goes to infinity. The curvature (i.e.

the reciprocal of the radius) goes to zero and the

isogonic curves become lines or rays which originate at

the detector. The exit paths of the scattered photons (or

neutrons) are, therefore, colinear with the isogonic

lines (see Appendix II). The same equations describing the

detector response for the wide beam source/uncollimated

detector case are valid provided the appropriate limits on

the integrals are used.








Discrete Forward Equation

This section deals with the practical aspects of

expressing equation (2.7) in a form which enables computing

the integral of equation (2.6*) using a digital computer.

The details of how the computations are done are presented

in Appendix III. The purpose here is to show the relation-

ship between the discrete equations, which are used in

reconstructing an image of a flow cross section, and the

"exact" integral form.

The measured data are available as a discrete set from

which the singly scattered flux spectrum must be extracted.

The extraction process of removing the multiply scattered

flux component and correcting for the detector system

efficiency can be expressed as



S1 = (D B )/c (2.8)



where S1 is the number of photons that arrived at the

detector having suffered only a single Compton scattering

collision. The notation Sl denotes a vector with components

Sc. The object here is to be able to calculate S1 given the

density distribution in the illuminated cross section of the

pipe. An expression for calculating Sc is derived using

equations (2.7) and (2.7*)



S = sOf dEf wPAiAodZ (2.9)
AEc M(E)









Beam
Collimators Pipe Wall



Monoenergetic
Source \
Array




Isogonic Line

Detector
Figure 2.6: A parallel beam illumination/uncollimated
detector gamma-ray scattering system.





Isogonic Segment
..i ....... ... iiiiiii

nI K




S.iiiiiiii.. .......... .......






2 8 99 .



1 2 3 n


Figure 2.7: An n x n pixel grid is overlayed the flow cross
section in order to facilitate calculation or
the singly scattered photon flux. Individual
pixels are designated 2by a single number k
ranging from 1 to K(=n ).








The first step in numerically evaluating equation (2.9)

is to overlay the illuminated flow cross section with an

n x n pixel grid as shown in Figure 2.7. Equation (2.9) is

rewritten as a sum of integrals over the individual pixels



S = so I dEf wPkAiAodk (2.10)
kc AEc Z(E)


where pk(x,y) is equal to p(x,y) for points inside pixel k

and is zero elsewhere. The notation kc indicates that the

sum is to be carried out only over those pixels in the

isogonic segment c. The integrals in equation (2.10) can

be expressed equivalently by



S = Sol fdxfdyPwPkAiAoPSFc(x,y) (2.11)
k


The term PSFc(x,y) is the characteristic function of the

isogonic segment c defined as



PSF (x,y) ={ 1 if (x,y) is inside isogonic segment c
c 0 if (x,y) is outside isogonic segment c.



Any real imaging system imposes a lower limit on the

spatial resolution below which variations in the density can

not be seen. If the size of the pixels are about the same

as the spatial resolution limit, then it is reasonable (even

necessary) to assume a constant density within each pixel.








The remaining integral in equation (2.11) is approximated by



S = S I pkAckPSFk (2.12)
k



The factor Ack is an average quantity defined by



ffdxdypw(x,y)Ai(x,y)A (x,y)
Ack = -k (2.13)
-k ffdxdy
k



where the integrals are carried out only over the area of

pixel k. The function PSFck is the fractional area of pixel

k which is inside isogonic segment c or



PSFck = ffdxdy PSFc(x,y). (2.14)
k



An implicit assumption in equation (2.12) is that the

attenuation factors do not vary significantly across a

pixel. This assumption is valid except for pixels near an

edge in the density distribution and simultaneously on the

edge of an isogonic segment. The assumption also breaks

down when the size of the pixels are on the order of a

photon mean-free-path. In practice there are no sharp edges

to the isogonic segments due to the scattering angle blur.

The scattering angle blur causes the isogonic segments to

"feather" off at the edges (and overlap) as shown in








Figure 2.8. The effect is that the assumption works in real

applications where non-point source and detectors are used.

The definitions for Ack and PSFck must be modified to

account for the scattering angle blur. The new definition

for Ack specifies that the attenuation factors must not only

be averaged over the area of the pixel but also over the

range of possible energies with which the photons can reach

the detector. The new definition of PSFck is the fraction

of photons reaching the detector from pixel k with

energies in the range of channel c multiplied by the pixel

area.

Equation (2.12) gives an expression for calculating the

singly scattered photon spectrum for a detector looking at

the scattered photon flux of a known density distribution

illuminated by a wide beam mono-energetic gamma-ray source.

Each channel of the singly scattered photon spectrum is seen

to represent an integral measure of the density distribution

within the corresponding isogonic segment. The goal of the

remainder of this report is to do the inverse operation of

extracting information about an unknown density distribution

given measurements of the singly scattered photon flux.

Chapter 3 tackles this inverse problem by correlating the

size and shape of spectra with the different possible flow

regimes. Chapter 4 goes one step further in solving the

inverse problem and looks at ways of inverting

equation (2.12).

























PSFch(x,y)


Figure 2.8: The scattering angle blur causes the isogonic
segments to feather off and overlap. PSFch(x,y)
is the fraction of those photons which scatter
at point (x,y) toward the detector that have an
energy in the spectrum channel ch.















CHAPTER 3
METHOD OF SPECTRAL MOMENTS



The basic premise of Compton scatter tomography is that

the singly scattered photon spectrum contains information

about the density distribution. If this is true then each

flow pattern should impart some unique characteristics to

the spectrum. Examination of the first few moments of the

spectrum is a simple pattern recognition technique employed

in this chapter to extract these unique characteristics.


Description of the Method

Consider a Compton scattering setup using a wide beam

source and uncollimated detector as illustrated in

Figure 3.1. Mono-energetic gamma rays leave the source and

undergo Compton scattering interactions inside the pipe.

Some of the scattered gammas are received by the detector

and the energy spectrum of the detected gamma rays is accu-

mulated in an MCA. The total count in each MCA channel is

proportional to the density of the material in the isogonic

segment for the channel. Note that the source beam is

stopped down slightly in order to reduce the amount of the

pipe wall which is illuminated.

The typical two-phase fluid problem consists of a fluid

and its vapor. The vapor is, in general, considerably less

























































Figure 3.1:


Compton scatter densitometry configuration.
Scattered radiation from the illumination beam
is received by the detector. (Reproduced by
permission. Copyright o Instrument Society of
America (ISA) 1985 [21].)








dense than the liquid. Therefore, an MCA channel represents

an integral measure of the liquid in the isogonic segment.

The singly scattered photon spectrum is a projection of the

liquid distribution onto the energy axis. The energy axis

can, by suitable transformation, be related to any spatial

axis intersecting the isogonic lines.

A convenient reference axis perpendicular to the

isogonic lines is the Y-axis. The transformation relating

the energy axis to the Y-axis is available from

equation (2.5)



Y(E) = r(9(E)) + yc(e (E)) (3.1)



Note that each detector has a characteristic Y-axis.

The fluid distribution in the pipe can be inferred from

the one dimensional image (i.e. projection) by examining the

moments of the spectrum with respect to the Y-axis. The

zeroth moment is a measure of the amount of liquid phase

present. The first moment is a measure of where the center-

of-mass of the liquid phase is. The standard deviation is a

measure of how the liquid phase region is distributed about

the center-of-mass. The moments are not direct indicators

of their respective quantities because of the attenuation

factors.

The zeroth moment is calculated in a straight forward

manner as










= Isc (3.2)
c



where Sc is a channel of the singly scattered photon

spectrum [see equation (2.8)]. The first moment and

standard deviation are calculated according to the formulae



Sc Y YdY

1 C AYc c+l
= cY Yc+l (3.3a)

1C
c



and



S C Y c (Y )2dY

2 c AY Yc+l
= (3.3b)

c sc
c



Figure 3.2 graphically depicts the zeroth moment, first

moment, and standard deviation for annular and inverted

annular flow at 50% void fraction. The two peaks in the

spectrum for the annular flow correspond to the liquid in

the illuminated portions of the annulus separated by the low

density vapor core region. The peak heights are not the

same because the source beam is attenuated before it reaches

the far side of the annulus which results in fewer scat-

tering collisions taking place there. Geometric attenuation





















1.00
z =0.313
M rel. units
0
0.75
0

0.50-

U T

0.25-



0 10 15 20
Y (cm)

(a)










Figure 3.2: Spectral moments of the singly scattered photon
spectrum. a) Spectrum for an annular flow;
b) Spectrum for an inverted annular flow.
(Reproduced by permission. Copyright ISA
1985 [21].)









































Y (cm)


(b)


Figure 3.2--continued









of the photon fluence (i.e. photons/area) for points further

from the source also occurs. Similarly, these same attenua-

tion terms are applicable with respect to the detector.

Only a single peak for the liquid core region is seen in the

inverted annular spectrum. The peak is not symmetric due to

attenuation as the gamma rays pass through the liquid and to

geometric attenuation.

The zeroth moments are not identical even though the

void fraction is the same. All of the liquid core region is

illuminated in the inverted annular flow. Not all of the

liquid annulus in the annular flow is illuminated. In

addition, the different liquid distribution leads to differ-

ences in attenuation, solid angle, and differential

scattering cross section. Consequently, the zeroth moment

of the inverted annular flow is larger than that of the

annular flow.

Both the annular and inverted annular flows are

symmetrically distributed within the pipe and so, the first

moment would be expected to lie at the center of the pipe.

However, due to the attenuation effects the first moment is

skewed closer to the source and detector. The standard

deviation of the annular spectrum is larger than that of the

inverted annular. This is expected since the annular

spectrum consists of two widely separated peaks while the

inverted annular spectrum is a more compact single peak.









Procedure and Results

Monte Carlo transport calculations are carried out to

test whether flow regime identification is possible by

examining the first few moments of the singly scattered

gamma ray spectrum. Several different flow regimes are

analyzed; in each case the energy spectrum from single

scattering interactions with liquid or vapor inside a pipe

is calculated for two detectors. The detectors are assumed

to be placed symmetrically with respect to the pipe and

source as shown in Figure 3.3. The source-pipe-detector

geometry is the same for all cases.

The source and detector are treated as points for

purposes of calculating the scattering angle. However, for

purposes of calculating the solid angle view of each

interaction, the detector has a radius of 0.5 cm. Detected

scattering interactions within the pipe wall are not

considered since the primary focus is to evaluate the gamma

rays coming from the fluid. However, attenuation due to the

pipe is taken into account.

The flow regimes considered are annular, inverted

annular, stratified, bubbly, and mist. The annular and

inverted annular flows are modeled as a liquid/vapor annulus

and a vapor/liquid central core concentric with the center

of the pipe, respectively. The stratified flow is modeled

as a pipe half-filled with liquid as illustrated in

Figure 3.4. The bubbly and mist flows are simulated by

choosing fifteen randomly sized (0.85








bubbles of vapor/liquid in a matrix of liquid/vapor, respec-

tively. New bubbles are periodically chosen during a cal-

culation to simulate stochastic flow conditions. The void

fraction is fixed at 50% for the annular, inverted annular,

and stratified flow. Due to the random nature of choosing

the bubbles, the bubbly flow has a void fraction of ~45% and

the mist flow has a void fraction of ~58%. Four calcula-

tions are made for each flow regime where the width of the

source beam in the plane of the pipe cross section varies

from 4.5 to 28 degrees. The results of the moment

calculations are tabulated in Table 3.1.

Interpreting the Moments

The moments are functions of the magnitude, shape, and

width of the spectra. The spectra change as the source beam

width varies since the amount of the flow cross section

illuminated changes. Stating general trends of the spectra

based on the changing amount of illuminated fluid is diffi-

cult due to the effects of attenuation. However, the

following two observations can be made. First, the magni-

tude of the spectrum from a narrow source beam is smaller

than that from a wider source beam. Assuming the source

strength remains constant, the number of photons illumina-

ting the pipe is less for the narrow source beam. The

result is that the zeroth moment becomes smaller as the

source beam narrows.

Second, the width of the spectrum decreases as the

source width is narrowed. As illustrated in Figure 3.5, the








Detector 1


Source '_,
ipe







Detector 2


Figure 3.3: One source, two symmetrically positioned
detector configuration used for the Monte Carlo
studies presented in this chapter. (Reproduced
by permission. Copyright ISA 1985 [21].)


Detector 1


Source


Vapor

Liquid


Figure 3.4:


4 Detector 2
Source-detector configuration relative to the
stratified flow regime case. (Reproduced by
permission. Copyright ISA 1985 [21].)












TABLE 3.1


Zeroth moment, first moment, and standard deviation of the singly scattered photon spectra
for the five flow regimes listed.
Source Beam Width


Flow Regime


ANNULAR
Detector 1
Detector 2

INV. ANNULAR
Detector 1
Detector 2

STRATIFIED
Detector 1
Detector 2

BUBBLY
Detector 1
Detector 2

MIST
Detector 1
Detector 2



[relative]


280

[cm]


0.311
0.313


0.337
0.340


0.621
1.00


0.434
0.436


0.431
0.423


<02>
[cm2]


140




0.0935 15
0.0935 15


3.7 0.194 16
3.7 0.193 16


6.1 0.360 17
5.4 0.353 15


8.4 0.204 16
8.4 0.205 16


6.6 0.183 16
6.9 0.184 16


80



15
15


3.7
3.6


0.0485 15
0.0484 15


0.108 16
0.108 16


4.50



3.5
3.6


6.6 0.207 16 6.9
6.3 0.158 16 6.6


7.6 0.108 16 8.1
7.5 0.107 16 8.0


6.5 0.108 16 5.9
6.7 0.112 16 5.7


0.0267 15
0.0267 15


0.0601 16
0.0600 16


0.116 16
0.0789 16


0.0597 16
0.0602 16


3.5
3.5


6.7
6.4


8.2
8.2


0.0614 16 5.3
0.0600 16 5.2


Note: The spectra are from a Cs-137 source illuminating a 12 cm diameter pipe with a 0.75
cm thick steel wall. The beam divergence normal to the flow cross section is 5
degrees in all cases. The void fraction is nominally 50% for all cases.









region within the isogonic segments associated with the

extreme channels of the spectrum is no longer illuminated as

the beam narrows. The effect is more pronounced for the

isogonic segments at low values of the Y-coordinate. If

liquid is present in the affected regions, the result is

that the spectrum narrows and the narrowing can be more

pronounced at the low end of the Y-axis. Such an asymmetric

change in the spectrum width means the first moment shifts

to larger values of Y as the beam narrows. A narrower

spectrum would also be expected to yield a smaller standard

deviation. The magnitude of change in the first moment or

standard deviation, or whether any change occurs at all,

depends on any accompanying change in the shape of the

spectrum. If there is no liquid in the affected regions

then, of course, there is no change in the width of the

spectrum.

A Detailed Example

The following is a discussion in physical terms of how

the moments vary with changes in the source beam width. So

that sufficient attention can be given to the changes in

each moment, the discussion is confined to only the strati-

fied flow regime. Enough insight should be gained to

interpret the data in Table 3.1 for the remaining flow

regimes.

Figures 3.6 and 3.7 are the spectra for the upper and

lower detector, respectively, viewing a stratified flow.

These figures show the progressive changes in the magnitude,









shape, and width of the singly scattered photon spectrum as

the source beam width narrows. The magnitudes of the spectra

decrease as the beam narrows because there are fewer photons

illuminating the flow. The magnitudes of the two figures

should not be compared since the vertical scales are diffe-

rent. Instead the figures should be compared on the basis

of the zeroth moments.

The zeroth moment of the lower detector's spectrum is

larger when the source beam is wide, while the zeroth moment

of the upper detector's spectrum is larger when the source

beam is narrow. The behavior of the zeroth moments can be

explained by considering the effect of attenuation on the

scattered gamma rays and the differences in solid angle each

detector views of the scattering interactions. When the

source beam is widest the average attenuation of the

scattered gamma rays reaching each detector is about the

same. The average position of scattering occurs in the

middle of the liquid and the scattered gamma rays must pass

through almost the same amount of liquid and steel to reach

either detector. When the source beam is narrowed the

average scattering interaction takes place toward the center

of the pipe and the scattered gamma rays must pass through

more liquid to reach the lower detector. Thus, attenuation

of the scattered gamma rays explains why the zeroth moment

of the upper detector is larger when the source beam is

narrow. When the source beam is wide the zeroth moments

should be the same since the average attenuation is the


























































Figure 3.5: Darkened regions are those isogonic segments
which are eclipsed when the more restricted
beam illumination is used. (Reproduced by
permission. Copyright O ISA 1985 [21].)






























Z beam line zeroth
3 width pattern moment
0 (degrees) (relative)
.75 28 0.621
0 22 ---- 0.532
'- 14 --- 0.360
S8 -- 0.207 \
1- 5 4.5 -- 0.116
LU


S.25-
h-
U.1 -- -


















scattered gamma rays for upper detector viewing
a stratified flow as in Figure 3.4. (Reproduced
by permission. Copyright ISA 1985 [21].)




























1.0
I-- beam line Izroth
2 width pattern moment
(degrees) (relative)
0
0 .7 28 1.000
S22 ---- 0.708
o 14 --- 0.353 A
-.. 0.158
0 4.5 ---- 0.079
F- .5- \ .
LU \


p.25-



0 5 10 15 20
Y (cm)















Figure 3.7: Detected photon energy spectra of once-
scattered gamma rays for the lower detector
viewing a stratified flow as in Figure 3.4.
(Reproduced by permission. Copyright ISA
1985 [21].)









same. The reason they are not is because the upper detector

is further away and presents a smaller solid angle to each

scattering interaction.

The behavior of the first moments can be explained in

the following manner. Observe that the spectrum for the

widest source beam is peaked at about 13 cm for the lower

detector and at about 17 cm for the upper detector. Due to

the source-pipe-detector geometry the isogonic segments do

not all intersect the liquid region. The projection of the

liquid region onto the corresponding Y-axis for the two

detectors is not symmetric. As the source beam narrows both

spectra become more uniform since approximately the same

amount of liquid is illuminated in each isogonic segment.

The remaining slight peak is due to the non-uniform amount

of attenuation and to solid angle considerations. The more

uniform a spectrum the more the center-of-mass coincides

with the center of the spectrum. The result is that the

first moment shifts to higher Y values for the lower

detector and shifts to lower Y values for the upper

detector. Note also the preferential narrowing of the

spectra which tends to enhance the changes in the first

moments.

The standard deviation for both detectors' spectra

increases uniformly as the source beam narrows from 28 to 8

degrees and then decrease slightly as the beam narrows

further to 4.5 degrees. When the source beam is widest both

spectra are peaked near the center-of-mass which gives a









relatively small standard deviation. The standard deviation

increases as the source beam narrows because the spectra

become more uniform. The slight narrowing of the spectra

leads to a smaller standard deviation. The two competing

effects result in the observed behavior.

Using the Moments

The purpose of calculating the moments is to identify

the flow regime automatically and without any a priori

knowledge. Inter-comparing the moments of all the flow

regimes at each particular source beam width suggests that

using the widest source beam offers the greatest discrimina-

tion between flow regimes. A close agreement between the

moments of the upper and lower detectors is seen in all the

data for the annular, inverted annular, bubbly and mist

flows where the fluid distribution is symmetric. The

stratified flow is asymmetric, but only for the wider source

beams (i.e. > 14 degrees) is the asymmetry discernible from

comparing the moments of the upper and lower detector. The

annular flow is discernible from the other symmetric flows

with the wide source illumination because of its large

standard deviation. The inverted annular flow conversely is

distinguishable because of its small standard deviation.

Discriminating between the bubbly and mist flows is not

possible on the basis of the first moment or standard

deviation due to their similar nature. Bubbly flow is

expected when the void fraction is low (i.e. < 50%) while

mist is expected when the void fraction is high. The zeroth









moment is sensitive to the void fraction and could be used

to determine which flow exists.

The results presented so far show that flow regime

identification can be made using the singly scattered photon

spectra from two detectors sensing the emerging scattered

photon flux from a flow cross section illuminated by a wide

beam source. A tomographic image of the flow can be

produced using a simple model of the flow cross section for

each flow regime. The consistency between the actual flow

regime and the flow regime as identified by the spectral

moments can be checked using the solution to the forward

problem [i.e. equation (2.12)]. The calculated singly

scattered photon spectra for the two detectors based on the

identified flow regime should match the measured spectra.

If the calculated and measured spectra are not close enough,

then the flow cross section model is adjusted accordingly.

For example, if the spectral moments indicate that the flow

is annular, then the size of the annulus is adjusted until

the measured and calculated spectra match. The next chapter

develops other methods of producing such tomographic images

of the flow. These other methods are more systematic in

adjusting the flow cross section image so that the measured

and calculated singly scattered photon spectra match.














CHAPTER 4
THE RECONSTRUCTION PROBLEM



A general statement of the problem at hand is


Given a set of singly scattered photon flux measure-
ments, {S }, made using a wide-beam source and
uncollimated detectors, reconstruct the liquid-vapor
phase distribution {p} that produced {S }.


Ideally the set of flux measurements is acquired with

a single exposure. That is, the source and detector

locations are fixed and there are no moving parts. Insight

into this problem can be gained by considering the following

formulation of the forward problem from equations (2.6*)

and (2.7):



SI=[T]p (4.1)



where [T] is the matrix which mathematically describes the

scattering transport process, S1 is the vector of singly

scattered photon flux measurements, and p is the vector of

pixel densities. The obvious solution of inverting [T],

that is,



p=[T]-1 S1 (4.2)








is not strictly possible since the elements of [T] are

functions of p. The problem at hand is nonlinear. Until

now the general approach has been to avoid dealing directly

with the nonlinearity of [T].

One often-proposed method (e.g. [24, 25]) is to use

multiple source-detector locations and/or source energies in

order to cancel out the nonlinear aspects of [T]. There are

numerous practical limitations in implementing this approach

for dynamic measurements. Another proposed method is to

assume a density distribution based on prior experience and

use this in computing [T]. Since [T] is a function of the

density, a bias in p is likely for the high contrast

reconstructions expected from two-phase flows. A charac-

teristic of these methods is the idea of calculating the

density distribution in an algebraic, once-through fashion.

It is proposed here that the problem be viewed as a

tomographic reconstruction problem where the measured data,

S, represent projections of the phase distribution, p, onto

non-physical energy axes. The projection data are

incomplete since only a few detectors are employed. Simila-

rly the elements of [T] are at best known only approximately

since they depend on the unknown phase distribution p. When

the available information is incomplete and/or imperfect, an

algebraic solution may not exist or several solutions may

exist. These are familiar situations in CT and the analytic

techniques are already developed to handle them. The analo-

gous projection operator, [T], of current CT applications is









mathematically well-behaved (i.e. linear, single-valued,

etc.). The CT reconstruction methods must, therefore, be

modified to handle the nonlinearities of Compton scatter

tomography (CST). Several CT reconstruction techniques are

adapted to CST in this chapter.


General Solution Methods


The basic methods of CT for solving the reconstruction

problem can be grouped into five categories:



1. Fourier transform: exponential Radon transform [26],
filtered backprojection [27].

2. Direct scanning [12, 14, 15].

3. Matrix inversion: successive approximation [20],
maximum entropy (see Appendix IV) [28], generalized
inverse, and pseudo inverse [29].

4. Series expansion [30]: algebraic reconstruction
technique (ART), simultaneous iterative reconstruc-
tion technique (SIRT), and iterative least squares
(ILS).

5. Statistical: Monte-Carlo backprojection.



The transform techniques are recognized as being the

more elegant, mathematical reconstruction methods. Trans-

form methods in general require a large number of projec-

tions and do not show good noise immunity. The cost of

detectors and the physical limitations in mounting and

cooling the detectors imply a limit on the number of projec-

tions available. Practical limits on source size, limits on

detector size, and desired framing rate are all going to








result in statistical noise in the measurements. These

features of using Compton scatter to image dynamic two-phase

systems make transform methods unsuitable as a recon-

struction technique.

Direct scanning techniques inherently require a narrow

collimated source beam and a moving source detector arrange-

ment. A narrow collimated source is less photon efficient

than a wide collimated source. This inefficiency coupled

with the time necessary to translate the source and detec-

tors means that a direct scanning system is slower than a

comparable wide-beam source, uncollimated detector system.

Since one goal is to measure dynamic two-phase flow, the

faster measurement system is the obvious choice. There is

another reason for dismissing direct scanning techniques

that involves accounting for the multiply scattered photon

flux. The structure of the flow regime is more strongly

imprinted on the multiply scattered photon flux when a

narrow-beam source is used than when a wide-beam source is

used. Any structure in the multiply scattered photon

spectrum means that approximating the multiply scattered

photon spectrum by a general function B(E) [see equation

(2.6*)] is not possible. Removing the multiple scatter

component in direct scanning systems, therefore, requires

more computation and is more prone to error.

Matrix inversion and series expansion reconstruction

techniques are equivalent in that a solution is sought by

iteratively improving an initial guess of p (and hence [T])









until some integral transformation of the reconstructed

phase distribution and the true phase distribution are

judged close enough with respect to the measurement error.

The difference between matrix methods and series expansion

methods is that at each iterative step matrix methods update

all the pixel values at once, whereas series expansion

methods update the pixel values sequentially. Matrix

methods usually exhibit better convergence rates than series

expansion methods since they can take into account higher

order effects in a more straightforward manner. Techniques

from either of these categories appear to be well suited for

use in a CST system.


A Closer Look

This section details how specific CT techniques are

adapted for use in CST. The first technique to be

considered is backprojection. Generally the initial step in

a reconstruction algorithm is a backprojection operation

which yields a very smooth representation of the true image.

Much of the fundamental differences between traditional CT

and CST are evident in the way the projection data are

handled in order to perform a basic backprojection opera-

tion. The second technique considered is Monte-Carlo

backprojection (an original concept by this author and Dr.

E. E. Carroll at University of Florida). The result of

Monte-Carlo backprojection is a higher contrast image than

obtained with traditional backprojection.









The third part of this section is a general discussion

of matrix inversion techniques. The purpose of this

discussion is to introduce the basic iterative nature of any

practical solution strategy imposed by the nonlinearity of

CST. The fourth part of this section adapts three series

expansion techniques (ART, SIRT, and ILS) for use in CST.

Backprojection

Consider a projection formed by the conventional CT

scanner shown in Figure 4.1. The pencil beam source irra-

diates the object and the transmitted radiation is measured

at points in the image plane. The source and detector are

moved in the direction indicated to scan the entire object

and form a single projection. Other projections are formed

by rotating the source and detector around the object and

taking another scan. The number of photons reaching the

detector at each point in a projection is a function of the

average optical path through the object. The functional

relationship is



Ijc = Ioexp{-f p(s)pt(s)ds} (4.3)
kjc


where Ijc is the output of the detector at position c in
projection j,

I is the number of photons leaving the source,

p(s) is the density at point s,

pt(s) is the total mass attenuation coefficient of the
material at s, and

Pjc indicates that the integral is to be carried out













Pencil Beam Source


















-- ----Image Plane


Detector






















Figure 4.1: Example of a conventional transmission CT
scanner. The source and detector move parallel
in the direction indicated by the straight
arrows, and also rotate about the object being
imaged.









along the line between the source and detector
when at position c in projection j.



In practice one does not work with Ijc, but instead with Rjc

given by



K
Rjc = -n(Ijc/Io) = t Pksjck (4.4)
k=l



where pk is the density in pixel k,

Asjck is the length of line jc contained in pixel k,
and

K is the total number of pixels.


The value of Asjck is nonzero only if line jc crosses pixel

k and so the sum in equation (4.4) is equivalent to the

integral in equation (4.3). Since the liquid and vapor are

of the same material, the mass attenuation coefficient is a

constant and is taken outside the integral. The average

density, Pjc, along line Zjc can be computed from equation

(4.4) as



K
Pjc = Rjc/{t Asjck} = Rjc/Ljc (4.5)
k=l



A backprojection image is made by forming a composite of the

average line densities, Pjc'









J C(j)
bpk ={ I PjcASjck/
j=l c=l


J C(j)
1 sjck
j=l c=l


(4.6)


where bpk is the density of pixel k in the reconstruction,

and Asjck is used as a weighting factor.

An analogous relationship for computing a backpro-

jection image in CST can also be developed. The quantity

corresponding to Rjc in CST is the singly scattered photon

flux values, Sc, given by


K
S= S PkAjkPSFjckAs
k=l


where p is the incoherent scattering coefficient,

S is the uninteracted source flux,

Pk is the density in pixel k,
As is the length of one side of a (square) pixel,

Ajk is an average quantity which takes into account
all aspects of photon transport from the source to
pixel k, and from pixel k to detector j, and

PSFjck is the.fraction of the photons arriving at
detector j from pixel k which have an energy in
the range of channel c (Ec to Ec+l).


Note that PSFjck is zero if it is not physically possible

for a singly scattered photon from pixel k to arrive at

detector j with an energy in the range Ec to Ec+1. So, the

sum of equation (4.8) represents an integral over the

isogonic segment c of detector or projection j and is

equivalent to the sum of equation (4.5). Another point to


(4.7)










note about equation (4.7) is the that double subscript jc

can be replaced by an equivalent single subscript g,



K
S = o PkAgkPSFgkAS (4.7*)
k=l

which is convenient in formulating the problem in matrix

notation as in equation (4.1).

The average density in an isogonic segment, Pjc, is

computed as



K
jc = Sc /AsS I AjkPSFjck (4.8)
k=l



and similarly the backprojection image densities, bpk, for

CST are computed as



J C(j)
I I PjcPSFjck
bpk = j=l c=l .(4.9)
J C(j)
1 l PSFjck
j=l c=l


A basic problem with equation (4.8) is that the density

distribution must be known in order to calculate the Ajk

factors. A reasonable solution is to use a homogeneous

distribution with an average void fraction that is either

assumed or gained from an independent measurement. The

average void fraction of a CST backprojection image must









then be checked for consistency with the initial void frac-

tion. Any inconsistency can be resolved by some iterative

refinement. Equations (4.5) and (4.6), in contrast, require

no knowledge of the average void fraction and a one-step

backprojection computation is possible in CT.

Monte-Carlo Backprojection

The motivation for using a statistical approach to form

a backprojection image comes from three observations.

First, the singly scattered photon flux, Sjc, is directly

proportional to the average density in the isogonic segment

jc. Second, any detected photons are due almost entirely to

scattering in the liquid phase since the density ratio of

liquid-to-vapor phase for water is usually 1000 to 1 in

nuclear reactors. Third, the features of one detected

spectrum can be cross correlated to the features of another

spectrum detected at another location. These observations

are summarized by the following example. If a single peak

is observed in channel c of detector j and a similar peak is

observed in channel c" of detector j', then the liquid phase

is localized to the region where isogonic segment jc

intersects segment j'c'.

Monte-Carlo backprojection treats each detected

spectrum as an independent one dimensional probability

density function of where the liquid phase is located. A

reconstruction of the liquid phase distribution is made by

statistically cross correlating the detector spectra. The

difference between Monte-Carlo backprojection and









traditional backprojection is that traditional backprojec-

tion ignores the cross correlation information.

The initial step in Monte-Carlo backprojection is to

form an independent probability density function (pdf) from

the modified average isogonic segment densities of each

detector

PiC
pdf (c) =. (4.10)
C(j)
SPjc
c=l



The average isogonic segment density is modified in the

following manner:



K
P;c = S/ Asso X Ajk (4.8*)
k=l



in order to account for the non-uniform spatial width of the

isogonic segments.

Assume that a fine-mesh square pixel grid is employed

and that a pixel is either liquid filled or vapor filled.

Provided a fine enough mesh is used, most density distribu-

tions can be reasonably represented by such a scheme. The

number of pixels, N, which are liquid filled is easily

computed from an estimate of the average void fraction, a,


N = round((l ()-K)


(4.11)









where K is the total number of pixels. All the pixels are

assumed vapor filled to start and the determination of which

N pixels to designate as fluid filled proceeds as follows:


1. Given a set of J random numbers {r-}, an associated
set of channel numbers {chj} is generated by
inverting the pdf 's

ch
rj 1 pdfj(c); j=l, 2, ..., J (4.12)
c=l


There are J detectors and one channel from each
detected spectrum is chosen by equation (4.12).


2. Associated with each channel chosen is a set of
pixels, Set(ch.), which contains all the pixels
within the corresponding isogonic segment. A new
set is formed by taking the intersection of all J of
these Set(chg)'s. The intersection set, I-Set,
contains a t6tal of M pixels. Taking a new random
number, rJ+1, the mth element


m = integer(M-rj+1) + 1 (4.14)


of I-Set is chosen as the pixel to be fluid filled.
If M is zero, then no pixel is designated as being
fluid filled.


3. Steps 1 and 2 are repeated sequentially until the
necessary N pixels are turned on.


Figure 4.2 compares the results of Monte-Carlo and

traditional backprojection. The Monte-Carlo reconstruction

is not as blurred and better represents the actual flow

regime. Although the reconstructed image resembles the

correct flow regime, a comparison of the calculated and

detected singly scattered photon flux spectra shows that

there is still room for improvement (Figure 4.3).
















(a)














(b)















(c)


META









ULFd


.Detector


Source *


-.Detector


Figure 4.2: Reconstruction of a 50% void fraction annular
flow regime. a) Monte-Carlo backprojection;
b) Simple backprojection; c) source and detector
locations.









Calculated -------
"Measured"


(a)


(b)


Figure 4.3: Comparison of the calculated singly scattered
photon flux for the Monte-Carlo backprojection
reconstruction of Figure 4.2 to the "measured"
spectra. a) Upper detector; b) Lower detector.










The addition of a feedback mechanism to the pixel

selection process is one possible improvement. Instead of

randomly choosing a pixel from the intersection set, those

pixels which improve the calculated spectra could be chosen

preferentially. A more complete evaluation of the Monte-

Carlo reconstruction method and possible variations is left

for future work.

Matrix Methods

Generally matrix inversion methods are avoided in CT

due to the large size (i.e. > 1000 x 1000) of the matrices

involved. Even in CST where there are only a few

projections and, with a coarser pixel mesh, the matrices can

be on the order of 200 x 200. Inverting matrices of this

size requires large amounts of storage, long computing

times, and is prone to round off error problems. Whatever

the difficulties in implementing a matrix method solution,

the compact notation is useful.

There are many ways of formulating the reconstruction

problem as a series of matrix operations. The most basic

approach is to treat the problem as an optimization problem.

A search is made for the density distribution, p*, which

minimizes, either directly or indirectly, the difference

between the measured data, S1, and the calculated singly

scattered photon flux, F, that is,

t All data in this report, unless otherwise stated, are
the result of Monte-Carlo photon transport calculations
(i.e. numerical experiments) and not physical
measurements.









min IS1 FI or min 1S1 [T(p)]p| (4.15a)
P P


subject to the constraints



Pvapor < Pk < Pfluid k=l,2,...,K. (4.15b)


The notation I..I represents some norm functional.

Any reconstruction algorithm in CST has a structure

similar to that depicted in Figure 4.4. Note that the

nonlinearity of the problem requires the recalculation of

the transport matrix for each new density distribution, pm.

Changes in the density distribution directly affect the

attenuation of the source photons and scattered photons

which is accounted for in the transport matrix. If the

variation in pm is small so that [T] is approximately

constant, then it may be possible to skip the recalculation

step for a few iterations. In contrast, the corresponding

matrix in CT is only a function of the system geometry and

remains constant throughout the calculation.

Another point at which the nonlinearity of CST needs to

be addressed is in the calculation of Em+l. One method of

calculating Tpm+l is to use the method of steepest descent


pm+l = (Jacobian of P(pm))-1.(l p( m)) (4.16)



The transport matrix can be treated as a constant or as a

function of the density in the calculation of the Jacobian.

























































Figure 4.4: Flow chart outlining the general structure of a
reconstruction algorithm for Compton scatter
tomography. The superscript m denotes the
iteration index.










Treating [T] as independent of p means a smaller step size,

f, is necessary and that Tp does not point directly toward

the path of steepest descent. Including the dependence

of [T] on p results in a significant increase in the storage

requirements and in computing time.

Algebraic Reconstruction Technique

ART is historically the first of the series expansion

techniques applied to commercial CT devices [31]. There are

two basic forms of ART algorithms: multiplicative and

addative. Presented here is a variant of the additive

algorithm.

A single ART iterative step involves correcting the

density vector by adding to it a correction vector



5m+l = pm + cm+l.wm+l (4.17)



where cm+l is a constant correction factor and wm+l is a

weighting vector. The correction factor is computed so that

one of the calculated singly scattered photon flux values is

consistent with the corresponding measured value, that is,



S Fm+l (= Tk p ) (4.18)




Subsequent iterations correct the density vector with

respect to the other measured values, one at a time. The

subscript g is incremented cyclically








gm+l = (gm + 1) modula G + 1



with each iteration in order to cover all G data points.

The correction factor can be derived by substituting

equation (4.17) into (4.18)

S1 VTm+Ipm
cm+l 9 g k gm+l. (4.19)

yTm+lwm+l
gk



The transport matrix elements are unknown until the new

density vector is known which in turn requires the

correction factor. One solution to this predicament is to

use the transport matrix calculated during the previous

iteration, that is,


ql TM Mpi
cm+l S gk P gm+l (4.19*)

ITm m+l
TgkWk



This approach is valid whenever the transport matrix does

not change much from one iteration to the next, as when

close to the final solution.

There are several possible choices of weighting

vectors. The simplest choice is an identity vector,

however, this results in all the pixel values being altered

including pixels which do not affect the sum in equation

(4.18). An identity vector modified such that









m+l 1 if Tg 0m+
w, =0 gk f m+l
k { 0 if Tgk 0 g = g (4.20a)



sensibly corrects only those pixels contained in the

isogonic segment corresponding to S1. Another reasonable

weighting vector is



w = Tk g = gm+1 (4.20b)



which effectively weights a pixel according to the proba-

bility that a scattered photon from that pixel is counted in

S A third possible weighting vector is



wm+1 = PSFgk; g=gm+l (4.20c)



which is analogous to the fractional ray length weighting

scheme of CT [32]. The weighting scheme chosen here and

elsewhere in this report is that of equation (4.20b).

Figure 4.5 graphically shows how an ART solution

proceeds for a two data point, two pixel situation. The

lines Lm and Lm correspond to the two simultaneous equations

given by



Lm S1 = Tml + Tm22 ; g = 1,2. (4.21)



The lines change as the density vector is corrected since

the transport matrix is a function of the current density

distribution iterate. The lines tend to stabilize as the











1.0

L
L= L L,

1 2 2 3





p2
04Li L- - L \

0.8 L 3
22 L1




0.6 0 -



P 1

0.4 -







pp
0.2






0 0.2 0.4 0.6 0.8 1.0

P,
Figure 4.5: Graphic demonstration of how an ART solution proceeds for a two data point,
{Li}, two pixel, {pi}, solution. Superscripts identify the iteration number.









final solution is approached, indicating that the transport

matrix is only varying slightly.

Simultaneous Iterative Reconstruction Technique

SIRT is similar to ART except that the density vector

is corrected using all the projection data, 91, simulta-

neously. SIRT attempts to combine a complete cycle of ART

iterations (i.e. G iterations) in calculating the correction

vector. The correction vector elements are given by



I (Tmk/a ){(S Fm)/( ~Tk
p+ = g k (4.22)
I (Tmk/og)
g


where a2 is the measurement variance. The correction vector
g
resembles a backprojection image of the error vector



m = S1 F" (4.23)



which can be seen by comparing equations (4.22) and (4.9).

Reference 30 indicates that some method of damping may

be necessary to avoid a solution which oscillates and fails

to converge. Once the correction vector is obtained it is

multiplied by a constant damping factor, f. The damping

factor is calculated so that the variance weighted euclidean

norm of the error vector is minimized



min S1 [Tm](pm + fm+l1m+l) (4.24)
f










The solution to equation (4.24) is


+1 (em/2 )(Tm Am+1 H
{m+l g/g)(gkk)}
f = _g k (4.25)

I{(Tm m+l }2
{ k TkAPk )/Gg}
g k


The proposed SIRT algorithm for CST is


Pm+l = mid {Pvapor' Pk + fm+lApm+' Pfluid} (4.26)

where the correction vector is given by equation (4.22) and

the damping factor by equation (4.25).

Iterative Least-Squares Technique

The ILS technique computes a correction vector which

minimizes the variance weighted euclidean norm of the error

vector. This can be expressed as


U = I{(Sl Fm TkAp+)/a }. (4.27)
g K


The correction vector which minimizes U is found by setting

the derivative of U with respect to the correction vector

equal to zero. This leads to the set of equations



pm+l(TmkT /a) = iTm-em/aG; k'=l,...,K (4.28)
k g g

Solving the system of equations of (4.28) would violate the

principal motivation of going to an iterative approach,








which is introduced to avoid the computational overhead and

pitfalls of matrix methods. The usual course of action is

to assume that all the the correction vector elements are

zero except the element corresponding to the current pixel

k'. Equation (4.28) becomes

Vrm em/a2
Apk+1 = g k k'=l,...,K (4.28*)

(ITmk/og)2
g


Equation (4.28*) overestimates the magnitude of the indivi-

dual elements of the correction vector which necessitates

the introduction of a damping factor as is done for SIRT.

The form of the ILS damping factor, which is calculated

according to equation (4.24) also, is


1{(em/a2)(Tm -pm+l
m+l {(eg/ g)(TgkPk)}
fm+l = -g k g. (4.29)

Tm Apm+l)/7 }2
g k


The proposed ILS algorithm for CST is

p+1 = mid vaporo' m + fm+1lm+ Pfluid} (4.30)


where the correction vector is given by equation (4.28*) and

the damping factor by equation (4.29).








Comparison of Reconstruction Techniques


The ART, SIRT, and ILS reconstruction techniques are

compared here in order to select one of them to use in a

detailed evaluation of CST. Each of the techniques are used

to reconstruct the two flow regimes shown in Figure 4.6.

The flow regimes are analogous to the "hot spot" and "cold

spot" phantoms used in CT studies. The criterion used to

judge the best technique is the euclidean distance between

the true density vector and the reconstructed density

vector. The rate of convergence is considered as a

secondary criterion.

The projection data (spectra) used for each technique

are generated by Monte-Carlo transport calculations. The

singly scattered photon spectra calculated for the two

detectors are shown in Figures 4.7 and 4.8. Table 4.1 lists

the pertinent "experimental" conditions.

The reconstructions of the hot-spot and cold-spot flow

regimes are shown in Figures 4.9 and 4.10, respectively.

All the reconstructions exhibit undesirable streaking

artifacts. The artifacts are due to the lack of edge infor-

mation from using only two views. The reconstructions of

the hot-spot are better than those of the cold spot. Intui-

tively this is expected by observing that the signal-to-

noise ratio is higher for the hot-spot spectra. The ART

results have more pronounced streaking artifacts than either

the ILS or SIRT results. The poorer performance of ART is

in keeping with similar CT phantom study tests done for













Detector


Source




Collimator







(dimensions in centimeters)


K5


*- 17
of -- 17


Detector


-.j x



=6
15


.j


Figure 4.6: Source and detector configuration used for the hot-spot and cold spot tests.











































Energy [keV]
(a)


400


400


Energy [keV]
(b)


Figure 4.7: The measured singly scattered photon spectra for the hot-spot test.
a) Spectrum from the upper detector; b) Spectrum from the lower detector.




















\ iI'

4J _P
r i1



U)

o 0
U U







0 I 0 0
0 400 0 400
Energy [keV] Energy [keV]
(a) (b)
Figure 4.8: The measured singly scattered photon spectra for the cold-spot test. The
ordinate scale is the same as for Figure 4.7. a) Spectrum from the upper
detector; b) Spectrum from the lower detector.


















TABLE 4.1


Data pertaining to the hot-spot and cold-spot experiment.



Detectors: Number of detectors 2
Detector type planar
Diameter 1 cm
Active thickness 0.1 cm
MCA channel width 1 keV


Source: Isotope Cs-137
Physical shape cylinder
Diameter 0.31 cm
Height 0.63 cm


Material Density:


Hot-spot
bubble density
bulk density

Cold-spot
bubble density
bulk density


1.155 g/cc
0.0 g/cc


0.0 g/cc
1.155 g/cc










Y


-x


-.1
* 1.


ILS


t.2


SIRT
















ART


Figure 4.9:


C";:


:'1



4


"- 7


Reconstructions of the hot-spot flow pattern
by the ILS, SIRT, and ART algorithms.









Y



x


ILS
















SIRT
















ART


Figure 4.10: Reconstructions of the cold-spot flow pattern
by the ILS, SIRT, and ART algorithms.









medical applications. There is no discernible difference

between the ILS or SIRT results.

The rate of convergence and accuracy of the reconstruc-

tion methods is shown in Figure 4.11. The discrepancy

function is defined as the euclidean distance between the

true flow regime and the reconstructed flow regime. The

spectra difference function is defined as the euclidean

distance between the measured spectra and the spectra calcu-

lated using the reconstructed flow regime. Both the discre-

pancy and spectra difference functions are normalized by

their values for the initial density distribution, which is

an all vapor condition. The convergence of all the recon-

struction algorithms for the cold spot test is similar. The

ILS and SIRT methods show better convergence properties than

the ART method for the hot spot test as judged by the

discrepancy function. The ILS and SIRT method are more

desirable than the ART method due to their consistent

performance. There is no significant difference in the

performance of the ILS or SIRT methods based on examining

the discrepancy or spectra difference functions. The choice

of whether to use either the ILS method or the SIRT method

appears to be one of personal preference.

This chapter successfully demonstrates that CT recon-

struction algorithms can be used with CST. As vital as the

reconstruction algorithms are, they are not sufficient for

implementing a CST system. The next chapter considers the

effect of noise in the spectra on the reconstruction process












































Iteration


* ILS
o SIRT
A ART


5
Iteration


(a)
Figure 4.11: Comparison of the accuracy and convergence of the ILS, SIRT, and ART
algorithms, a) Hot-spot flow regime; b) Cold-spot flow regime.





























cu
a)










0
0)






0-






0)
4-
4-1





**-


Iteration


OILS
o SIRT
A ART


5 10
Iteration


Figure 4.11--continued.


>1
u


)
a4

0
Ul
-H















0




Full Text


TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT
USING COMPTON SCATTERING OF GAMMA RAYS
By
DAVID E. BODETTE
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1986

Copyright 1986
by
David E. Bodette

ACKNOWLEDGMENTS
I thank Dr. Alan M. Jacobs, my adviser, for his patient
guidance and careful attention to detail.
This work was partially supported by a Liquid Metal
Fast Breeder Reactor Technology Research Award from the
Argonne National Laboratory under contract with the United
States Department of Energy (DOE). The author was also
supported by the DOE Energy Graduate Traineeship Program.
iii

TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS iii
LIST OF TABLES vi
LIST OF FIGURES vii
ABSTRACT xi
CHAPTERS
1 INTRODUCTION 1
2 THE FORWARD PROBLEM 10
Collimated Point Source/Collimated Detector . . 11
Collimated Point Source/Uncollimated Detector 15
Uncollimated Point Source/Uncollimated Detector 17
Parallel Beam Source/Uncollimated Detector. . . 28
Discrete Forward Equations 29
3 METHOD OF SPECTRAL MOMENTS 3 5
Description of Method 35
Procedure and Results 42
Interpreting the Moments 44
A Detailed Example 46
Using the Moments 52
4 THE RECONSTRUCTION PROBLEM 54
General Solution Methods 56
A Closer Look 58
Backprojection 59
Monte-Carlo Backprojection 64
Matrix Methods 69
Algebraic Reconstruction Technique 72
Simultaneous Iterative
Reconstruction Technique .... 76
Iterative Least-Squares Technique 77
Comparison of Reconstruction Techniques .... 79
4 COMPTON SCATTER TOMOGRAPHY 90
An Example 90
Accuracy 92
Resolution 106
Image Contrast 109
IV

Page
6 SUMMARY AND CONCLUSIONS 118
Summary 119
A Practical Application 121
Conclusions and Future Research 126
APPENDICES
I MONTE-CARLO TRANSPORT PROGRAM 130
II INSIGHTS AND FINER POINTS 133
Parallel Versus Diverging Source Beams 133
Multiple Scatter Versus Detector Position . . . 137
Source Placement 142
Convergence Criterion 143
III RECONSTRUCTION COMPUTER PROGRAMS 144
Running the Programs 144
Flow of Programs 145
IV MAXIMUM ENTROPY METHOD 146
V DATA 149
REFERENCES 191
BIBLIOGRAPHY 195
BIOGRAPHICAL SKETCH 196
v

LIST OF TABLES
Page
2.1 Multiply scattered photon spectra 22
3.1 Spectral moment results 45
4.1 Data for hot-spot and cold-spot tests 83
5.1 Test parameters and values 98
5.2 Summary of reconstruction results 105
5.3 Parameters affecting image resolution 107
vi

LIST OF FIGURES
Page
2.1 Collimated source/collimated detector
configuration 12
2.2 Collimated source/uncollimated detector
configuration 16
2.3 Uncollimated source/uncollimated detector
configuration 19
2.4 Illustration of scattering angle blur 25
2.5 Illustration of an isogonic segment 27
2.6 Parallel beam source/uncollimated detector
configuration 30
2.7 Overlaying of pipe cross section with an n x n
pixel grid 30
2.8 Feathering of isogonic segments 34
3.1 Compton scatter densitometry configuration 36
3.2 Spectral moments of a single-scattered spectrum 39
3.3 Source-pipe-detector configuration 44
3.4 Illumination of stratified flow regime 44
3.5 Pipe illumination for various beam widths 48
3.6 Spectrum of upper detector versus beam width 49
3.7 Spectrum of lower detector versus beam width 50
4.1 Conventional CT scanner 60
4.2 Monte-Carlo backprojection 67
4.3 Spectra from Monte-Carlo backprojection 68
4.4 Flow chart of reconstruction algorithms 71
4.5 Example of an ART solution for two pixel image 75
4.6 Hot-spot and cold-spot test configuration 80
4.7 Hot-spot spectra 81
4.8 Cold-spot spectra 82
vii

Page
4.9 Hot-spot reconstruction 84
4.10 Cold-spot reconstruction 85
4.11 Convergence of reconstruction algorithms 87
5.1 Calculated spectrum using incorrect source
strength value 91
5.2 Reconstruction using incorrect source
strength value 93
5.3 Experimental setup 95
5.4 Flow regime models 99
5.5 Singly scattered photon spectrum 100
5.6 Reconstruction with and without noise 101
5.7 Spectra with noise added 103
5.8 Detected spectrum with pipe wall effects 117
6.1 Spectrum from a 4 inch pipe 122
6.2 Detail of spectrum from the pipe 123
6.3 Effect of signal-to-noise ratio 125
6.4 Detector response versus bubble size 127
1.1 Example geometry of photon trajectory 132
11.1 Detector response for diverging source beam 135
11.2 Detector response for parallel source beam 136
11.3 Detector configuration for multiple scatter
versus detector position 138
11.4 Spectra for multiple scatter versus detector
position 139
V.l Spectra from annular flow 151
V.2 Spectra from annular flow . 152
V.3 Spectra from annular flow 153
V.4 Spectra from inverted annular flow 154
viii

Page
V.5 Spectra from inverted annular flow 155
V.6 Spectra from inverted annular flow 156
V.7 Spectra from bubbly flow 157
V.8 Spectra from bubbly flow 158
V.9 Spectra from bubbly flow 159
V.10 Spectra from stratified symmetric flow 160
V.ll Spectra from stratified symmetric flow 161
V.12 Spectra from stratified symmetric flow 162
V.13 Spectra from stratified asymmetric flow 163
V.14 Spectra from stratified asymmetric flow 164
V.15 Spectra from stratified asymmetric flow 165
V.16 Annular flow pattern 166
V.17 Reconstruction of annular flow 167
V.18 Reconstruction of annular flow 168
V.19 Reconstruction of annular flow 169
V.20 Inverted annular flow pattern 170
V.21 Reconstruction of inverted annular flow 171
V.22 Bubbly flow pattern 172
V.23 Reconstruction of bubbly flow 173
V.24 Reconstruction of bubbly flow 174
V.25 Reconstruction of bubbly flow 175
V.26 Stratified symmetric flow pattern 176
V.27 Reconstruction of stratified symmetric flow 177
V.28 Stratified asymmetric flow pattern 178
V.29 Reconstruction stratified asymmetric flow 179
V.30 Example of smooth convergence 180
IX

Page
V.31 Example of divergence 181
V.32 Example of divergence 182
V.33 Example of oscillating divergence 183
V.34 Example of non-convergence 184
V.35 Setup for laboratory experiments 185
V.36 Experimental spectrum of annular flow 186
V.37 Experimental and Monte-Carlo spectra 188
V.38 Experimental spectrum of stratified flow 189
V.39 Experimental spectrum of bubbly flow 190
x

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT
USING COMPTON SCATTERING OF GAMMA RAYS
By
DAVID E. BODETTE
December, 1986
Chairman: Dr. Alan M. Jacobs
Major Department: Nuclear Engineering Sciences
Using Compton scattered gamma-rays to measure local
void fraction was first suggested by Kondic and Hahn in
1970. The Compton scattered gamma-ray densitometers they
suggested employ a single, narrow beam source and either a
well collimated detector or an uncollimated detector. The
collimated detector configuration only gives the local void
fraction measurement in the small volume where the source
and detector collimators intersect. The uncollimated
detector configuration is a more efficient design since the
local void fraction along the source beam's path is measured
in a single reading. A logical extension of the technique
is to use wide beam illumination and uncollimated detectors
to sample an even greater portion of the flow cross section.
This report investigates and demonstrates several
methods of inferring two-phase flow parameters using wide
beam illumination coupled with two detectors placed
xi

symmetrically about the source and pipe. The spatial
distribution of the fluid in a slice of the pipe is encoded
with respect to energy in the singly scattered photon flux.
Two basic techniques are detailed for decoding the spatial
information: the method of spectral moments and the method
of computed tomography (CT).
Examination of the low-order moments of the singly
scattered photon spectra from the two detectors provides
sufficient information for flow regime identification. Flow
asymmetries are revealed by comparison of the spectral
moments of the two measured spectra. Further classification
of the flow pattern is made on the basis of the first and
second moments of the spectra.
The real focus of this report is the adaptation and
successful demonstration of CT techniques with Compton
scattering. Three series expansion techniques are
considered: the algebraic reconstruction technique (ART),
the simultaneous iterative reconstruction technique (SIRT),
and the iterative least squares (ILS) technique.
Application of the modified ART, SIRT, and ILS algorithms to
a hot-spot and a cold-spot reconstruction problem indicates
that SIRT and ILS are more accurate than ART. A more exten¬
sive testing of the ILS algorithm for a variety of model
flow regimes demonstrates the potential of Compton scatter
tomography as a quantitative measurement technique.
xii

CHAPTER 1
INTRODUCTION
The case for developing improved non-intrusive measur¬
ing technologies for sensing material/phase distribution in
two-phase flow has been well argued by many authors. The
problem is one which spans a broad number of fields inclu¬
ding manufacturing, engineering, security, and medicine.
Each potential application imposes its own peculiar demands-
which are often at odds with constraints in other applica¬
tions. This situation has resulted in many diverse
techniques being pursued and often optimized for each
specific problem. Examples of the technologies employed
include ultrasonics, microwaves, lasers, optics, and irra¬
diation with particles, X rays, or gamma rays.
One specific problem encountered in the nuclear power
industry and its associated research industry is to measure
the void fraction and/or liquid distribution (flow regime)
of two-phase flows inside steel walled pipes. Any intrusive
instrument which disturbs the flow pattern is certainly not
desirable, but an even more strict requirement often encoun¬
tered in these circumstances is that penetrations through
the pipe wall are not permitted. Another requirement of any
measuring instrument is that its output should be unambi¬
guous. In the case of two-phase flow void fraction
1

2
measurement this means the entire flow cross section must be
interrogated. There may also be physical limitations in
terms of size and orientation, as well as access to the
pipe.
A perusal of the recent literature on measurements in
two-phase flow [1, 2, 3, 4] reveals that the choice of
techniques for making non-intrusive local void fraction
and/or phase distribution measurements is limited to irra¬
diation with particles, gamma rays, and laser light. When
intrusion is not an issue, a variety of probe techniques are
available for local property measurements. Reference [4]
lists four basic probe methods which appear to be reliable
enough for general use: electrical probes, optical probes,
thermal anemometers, and micro-thermocouples. Laser
scattering and interferometry is suggested as a very
promising method when the flow boundary is not opaque.
There are techniques for making average void fraction
measurements which might be adaptable to the local void
fraction problem. Ultrasonics is capable of measuring the
average void fraction in very thick walled pipes, but only
for low void fractions (i.e. < 0.20) [3], An array of
ultrasonic sensors around a pipe might produce an image of
the flow provided the voids are not too big or too many.
Microwave resonant cavities are also capable of measuring
average void fraction. Using microwaves to sense local void
fraction would probably run into difficulties from the metal
pipes [1]. Nuclear magnetic resonance (NMR) is currently

3
receiving a lot of attention in the medical field and
certainly has potential for two-phase flow. The medical NMR
machines involve a radio-frequency field and a large
magnetic field which cause problems in nearby electronic
equipment. The metal pipes seem to be another source of
potential problems in implementing NMR for local .
measurements; however, plastic pipes can be used to avoid
such problems [5],
Gamma-ray attenuation is the predominant method of
choice for making average void fraction measurements and is
often used as the calibration standard for other methods.
Gamma-ray attenuation instrumentation and techniques are
recognized as having reached the highest level of
development. The reason for this dominance is that the
method places a minimum of constraints on the user. The
method is non-intrusive, fairly insensitive to flow pattern,
and functions over a broad range of void fractions, flow
regimes, and temperatures.
There are many ways of utilizing gamma rays or neutrons
to measure local void fraction and phase distribution. The
basic methods of implementing tomographic two-phase flow
measurement with gamma rays or neutrons can be grouped into
the three following categories:
1) Measure the emission of gamma rays or neutrons from-
activation products in the fluid, induced
excitation, or injected radiotracers.
2) Measure the attenuation of narrow or wide source
beams in traversing the pipe.

4
3) Measure the singly scattered gamma rays or neutrons
from a primary source beam.
Methods from category 1 could be roughly classified as
emission computed tomography (ECT) techniques. An example
of a category 1 method is a single photon ECT system which
examines the ~6 MeV gamma rays from the activation-decay
scheme pair
016(n,p)N16
n16 + 016 + Y + 0“.
Although recent advances in the analytical techniques of ECT
and in instrumentation suggest this as a viable method, it
appears that only cursory attention has been given to the
technique [6]. An ECT system could also be employed using
injected radiotracers; however, the quantity of radio¬
nuclides required poses a serious drawback [7].
The most familiar tomographic instruments comprise
category 2 and are generally referred to as computed axial
tomography (CAT) scanners. The basic hardware setup of CAT
scanners involves placing a source on one side of the pipe
and a detector on the opposite side of the pipe. The whole
source-detector assembly rotates in a controlled manner
around the pipe while the attenuation of the source beam at
different orientations is recorded. The rate at which CAT
scanners can rotate about the pipe is a limiting factor.
CAT scanners of the newest generation do not employ any

5
moving parts, but instead use a magnetically focused flying
spot electron beam and a fixed circular detector array [8].
These fifth generation CAT scanners can acquire an image in
approximately 1/10 second and are very promising where cine
tomography is necessary. The capital investment in these
new CAT scanners is significant and combined with the fact
that they are not compact or portable means they are not a
viable alternative in many applications.
Category 3 techniques are closely related to ECT
except that the emitted radiation does not come from
radioactive material within the fluid, but from radiation
scattered out of an external illuminating source beam.
Gamma-ray and neutron sources are possible candidates, and
methods based on both are actively being developed [9].
The earliest reported efforts with scattered radiation
focused on gamma rays scattering and involved measuring the
local density in a small volume [10]. A well collimated
source illuminated an object and a well collimated detector
measured the Compton scattered gamma rays coming from the
small volume where the source and detector collimators'
fields-of-view intersected. The output of the detector was
proportional to the electron density in the small volume
which in turn is proportional to the material density for a
constant material composition. The Compton scattering
method of local density measurement developed, at first,
independently in the medical and industrial imaging fields.
Lale published two articles, one in 1959 [11] and another in

6
1968 [12], which demonstrated how internal organs could be
imaged by exploiting the Compton scattering of gamma rays.
In 1970, Kondic and Hahn [13] published the first suggested
engineering application of Compton scattering, that of
making local density measurements in two-phase flow. They
proposed both a collimated source/collimated detector (as
Lale) and a collimated source/uncollimated detector arrange¬
ment. The uncollimated detector offered the improvement
that local density values were given for points along the
illuminated path of the source beam from a single
measurement. Acting independently, Farmer and Collins [14,
15] proposed and demonstrated the collimated source/uncolli¬
mated detector arrangement in determining anatomical cross
sections. Their images were only qualitative since the
calibration procedures used failed to completely remove the
effects of multiple scattering or variations in attenuation.
The analytic method for making quantitative measurements
using this source-detector arrangement by deconvolving the
substantial multiple scatter blur was developed by Anghaie
[16, 17].
The use of Compton scattering to make local density
measurements in two-phase flows was pursued earnestly in the
late 70's and in to the 80's. Almost exclusively the focus
was on the collimated source and either collimated or uncol¬
limated detector arrangement. Tomographic images produced
with either of these arrangements required scanning the
source and/or detector to cover all the points in the plane

of interest. Kondic [18] in 1978 first proposed using an
uncollimated or wide beam source coupled with an uncolli¬
mated detector to image an entire flow cross section. Such
an arrangement made more efficient use of the source photons
than the previous source-detector arrangements and required
no moving parts. Both of these attributes made the arrange¬
ment more adapted to the high framing rates necessary for
transient analysis. However, no critical study or viable
implementation of the wide source beam/uncollimated detector
arrangement has appeared in the literature.
Imaging with scattered neutrons was developed in the
late 7 0's and early 8 0's independently of and somewhat concur¬
rently of the present work with scattered gamma rays. The
use of neutrons to make local void fraction measurements was
suggested by Banerjee and Lahey [19] and separately by
other researchers. Neutron scattering tomography was
developed by Hussein [20] as his Ph.D. thesis. He used a
wide parallel beam of 14 MeV neutrons to illuminate a
simulated two-phase flow and then measured the fluence and
energy of the scattered neutrons using small organic scinti-
lator detectors at different positions around the flow
volume.
The focus of the present dissertation is the develop¬
ment and demonstration of the algorithms for reconstructing
the local density distribution of a two-phase flow using
Compton scattered gamma rays. The specific experimental
setup considered uses a wide-beam collimated gamma-ray

8
source to illuminate the flow and energy sensitive detectors
to examine the emerging singly Compton scattered photon
flux. To date, the reconstruction techniques proposed for
this setup have been based on deterministic algebraic solu¬
tions and have not addressed the problems of multiple scat¬
tering, statistical error, image blur, and other errors.
The contribution of this author is to apply the powerful
reconstruction algorithms developed in the medical and
optical imaging fields to this problem and to critically
assay the effects of„multiple scattering and statistical
errors on the reconstruction.
Tomographic imaging with Compton scattered gamma rays
has two distinct features which distinguish it from
previously encountered tomographic reconstruction problems.
First, the imaging processes dealt with in the medical and
optic fields have been linear or could be linearized by an
r
analytic technique. Imaging the phase distribution in a
two-phase flow with Compton scattered gamma rays is a non¬
linear problem. Second, the combination of using a point
source and the Compton scattering process introduces a
unique curved geometry to the flow reconstruction process.
The forward problem of predicting the emerging singly
Compton scattered photon flux given the local density dis¬
tribution is discussed in Chapter 2. Chapter 3 details a
method of inferring tomographic information about a two-
phase flow by examining the spectral moments of the Compton
scattered photons [21]. Chapter 4 develops several more

sophisticated tomographic reconstruction techniques and
compares their performance under ideal conditions using data
generated by Monte-Carlo calculations. Chapter 5 evaluates
Compton scatter tomography under more realistic conditions
and discusses some of the considerations in designing a
system. Chapter 6 summarizes the significant results of this
study and offers directions for future work.

CHAPTER 2
THE FORWARD PROBLEM
Fundamental to the problem of reconstructing the phase
or density distribution using Compton scattered gamma rays
is being able to model the gamma ray interaction with the
density field. The aim of this chapter is to describe the
basic equations predicting the emerging scattered photon
flux from a known density distribution illuminated by a
point qamma-ray source. The basic concepts of using gamma-
ray scattering to make local density measurements are
introduced first with the collimated point source/collimated
detector arrangement. Once the equations for analyzing this
setup are developed, the analysis is extended to collimated
point source/uncollimated detector arrangements by introduc¬
ing the equivalent approach of energy collimation. The
analysis is extended further to the wide-beam collimated
point source/uncollimáted detector arrangement. Lastly,
parallel source beam illumination is considered and the
differences between it and the diverging wide-beam
illumination are discussed. Throughout, a detailed analysis
is only given to the singly scattered photon flux. The
presence of the multiply scattered photon flux is included
in all the equations and discussed in general terms.
However, for the case of the wide-beam collimated source
10

11
/uncollimated detector configuration a specific method is
given for removing the multiple scatter component.
Collimated Point Source/Collimated Detector
Figure 2.1 illustrates the arrangement of a collimated
source and collimated detector for making local density
measurements. Monoenergetic gamma rays leave the source and
are collimated to a thin pencil beam which illuminates a
chord of the pipe. A similar collimator on the detector is
aligned so that its field-of-view intersects the source
beam at the small volume denoted by V. The detector
senses the source gamma rays scattered in the volume V.
The detector response, D(E), is
D(E) = eS°puwAiA0Az + M(E) (2.1)
where S° is the source activity;
p is the density of the fluid in the volume V;
pw is the incoherent differential scattering cross
section;
Az is the path length of source photons in V;
and AQ are the total attenuation functions, including
interaction attenuation and the solid angles of
view, along the incoming and outgoing paths
respectively;
e is the detector efficiency; and
M(E) represents the contribution to the spectrum of
multiply scattered photons.

12
Figure 2.1: A collimated source/collimated detector gamma-
ray scattering configuration for measuring the
density in the volume V.

13
The functional form of both attenuation functions is
Ax = Cx•exp(-J pxp(z)dz) (2.1a)
where Cx contains the geometric or field-of-view factors,
Px is the mass attenuation coefficient,
p(z) is the local density at z, and
x is a dummy subscript representing i and o for the
incoming and outgoing paths, respectively.
The line integral in equation (2.1a) is carried out either
along the incoming path, 1^, from the source to the volume V
or along the outgoing path, lQt from V to the detector. The
value of the integral is in units of optical path length.
The predominant interaction of gamma rays with energies
between 0.1 and 1.0 MeV in water is Compton scattering. The
energy of the singly Compton scattered photons arriving at
the detector is related to the polar scattering angle, ©, by
the kinematic relation
Eq = E±/{1 + (Ei/mec2)(l - cos 0)) (2.2)
where EQ is the energy of the scattered photon,
Ej_ is the energy of the source photon, and
mec is the rest mass energy of an electron.
The singly scattered photons, represented by the first
expression on the right hand side (RHS) of equation (2.1),
all arrive at the detector after scattering through

14
approximately the same angle and thus have approximately the
same energy. Assuming the attenuation along the incoming
and outgoing paths remains constant in time, the detector
response to the singly scattered photons is proportional to
the density in the small volume V.
The multiply scattered photons arriving at the
detector, unlike the singly scattered photons, have a broad
range of energies since they can suffer any number of
collisions (with concomitant changes in energy) prior to
reaching the detector. The functional form of the multiply
scattered photon component of the detector response is
M(E) = ////$(r ',E ',Q ' )K(r ',E ',Q ' ;r^,E,Q)dr 'dE 'dfi 'dr^dft (2.1b)
where $(..)is the angular photon flux and K (. -) is the
transport kernel. The kernel can be viewed as the probabi-
\
lity that a photon at r' going in direction ft' with an
energy E' suffers a scattering interaction at r'; the scat¬
tered photon goes in direction ft toward the detector with an
energy E; the scattered photon reaches point r^ in the
detector's active volume, interacts with the detector and is
counted. The integral with respect to r' is carried out
only over the detector collimator's field-of-view. The
detector collimator serves not only to pick out the volume
of interest, V, but to reduce the magnitude of the multiply
scattered photon spectral component by restricting the

15
region of space where the photons can suffer their last
collision prior to reaching the detector.
Collimated Point Source/Uncollimated Detector
The function of the detector collimator in picking out
the small volume V can also be achieved by electronic colli-
mation. An example is to use a high energy resolution
detector coupled to a multi-channel analyzer (MCA) instead
of the detector collimator. The point of scattering of a
singly scattered photon is related to its energy by equa¬
tion (2.2) and the polar scattering angle as determined from
the measurement geometry. The singly scattered photons
arriving at the detector now have a range of energies. The
detector response becomes
D(E) = S°p(x) vi^íxjAQUJeíE) + M(E) (2.3)
where S°, p, yw, Aj_, AQ/e , and M are as defined for equa¬
tion (2.1) and x is the point of scattering. The x-axis is
colinear with the source beam illumination path (Figure
2.2). Note that x and 9 can be expressed in terms of E
using equation (2.2) as
x = x^ -
l - cos^e
cos 9=1+ (mecz/E^)(1 - E^/E)
(2.4)

16
Source
Pipe
Configuration of a collimated source/uncollimated
detector gamma-ray scattering densitometer.
Figure 2.2:

17
where and are the detector coordinates illustrated in
Figure 2.2. The first term on the RHS of equation (2.3) is
strictly a function of a single variable for a fixed
geometry.
The multiple scattering term in equation (2.3) has the
same functional form as equation (2.1) except that the inte¬
gration with respect to r' is over all space since there is
no detector collimator. Therefore, the magnitude of the
multiply scattered photon component of the spectrum is
greater with an uncollimated detector than with a collimated
detector. Offsetting this increased measurement structured
noise is the fact that the singly scattered spectrum now
contains the local density information for the entire source
beam illumination path. The information which would have
required many measurements with a collimated detector is
available from a single measurement with an uncollimated
detector.
Uncollimated Point Source/Uncollimated Detector
Just as eliminating the detector collimator increases
the information contained in a single measurement, so does
eliminating the narrow source collimator. An uncollimated
or a wide angle-of-view collimated source can illuminate the
entire flow cross section. The detector response then
contains the local density information from every point in
the flow cross section. The manner in which the local
density information is encoded in the detector response,

18
however, complicates the data processing necessary to
reconstruct the local density distribution.
Figure 2.3 illustrates a wide beam source Compton scat¬
tering densitometer setup. The curves drawn across the pipe
represent lines of constant scattering angle, called
isogonic lines by Kondic [18]. Primary source photons
scattering at points along an isogonic line all scatter
through the same angle in order to reach the detector.
Using elementary geometry it is easy to show that the
isogonic lines are all segments of circles which pass
through both the source and the detector. The axis labeled
Y in. Figure 2.3 is the loci of all the centers of these
circles. The radius, r, and the location of the center-
point, yc, of each isogonic circle as a function of the
scattering angle, 9, are
r = I SD sin ©
2
yc = - -2SD cot ©
(2.5)
where SD is the distance between the source and the
detector.
The detector response is now represented by an integral
equation
ME)
D ( E ) = /
D(x,y)d£ + M(E)
(2.6)

19
ELEVATION VIEW
Figure 2.3:
Configuration of an uncollimated source/
uncollimated detector gamma-ray scattering
system for tomography.

20
where M(E) is the multiply scattered photon spectrum,
D(x,y) is the contribution to the detected photon
spectrum of singly Compton scattered photons
scattered at point (x,y), and
d£ is the differential line element along the
isogonic contour £(E).
The first term on the RHS of equation (2.6) is essentially
the integral of the first term on the RHS of equation (2.3).
The second term, as before, represents the multiply
scattered photon contribution to the spectrum and has the
same functional form as equation (2.1b).
The multiple scatter component can comprise approxi¬
mately one-third of the total signal [22, 23] and must be
removed. It may be possible to remove the multiple scatter¬
ing component analytically since the transport kernel is a
known function and the angular photon flux can be generated
from an order-of-scattering calculation once a good estimate
of the density distribution is obtained. However, this
would be a cumbersome process and does not appear to be
necessary as there appears to be an acceptable alternative
approximation.
A multiply scattered photon arriving at the detector
with energy E could have come from any point within a large
volume of the pipe. The structure of the fluid distribution
imposed on the multiply scattered photon spectrum would tend
to be averaged out because of the many possible paths that
the photon could have taken prior to arriving at the
detector with energy E. The result is that for wide beam
source/uncollimated detector arrangements the multiply

scattered photon spectrum should be smooth and a simple
background subtraction procedure can be employed to remove
it.
The validity of this intuitive argument has been
tested by a series of Monte-Carlo calculations (see Appendix
I). The multiply scattered photon spectrums have been com¬
puted for five flow regimes at both a high and low void
fraction. The calculations assume that the flow is
contained in a pipe whose thickness and composition are such
that it can be neglected. The results, summarized
pictorially in Table 2.1, show that to a first order approx¬
imation equation (2.6) can be rewritten as
D (E) = f D (x, y ) d£ + B ( E) .
2- (E)
where B(E) is some function of energy which is fitted to
the smooth multiply scattered photon
spectrum, and
D(x,y) is the contribution to the detected photon
spectrum of singly Compton scattered photons
scattered at point (x,y).
The function D(x,y) can be written more explicitly as
D(x,y)- S°e(E)yw(Ei,E0)p(x,y)Ai(x,y)AQ(x,y)
(2.7)
where e(E) is the detector efficiency,
Mw is the incoherent differential scattering cross
section for water,
p(x,y) is the local fluid density at (x,y),

22
TABLE 2.1
Monte-Carlo multiple scatter calculation results.
Flow Regime Void Fraction =0.7
Stratified
â–  D
Bubbly
â– D
r
Jy l
Annular
â– d
Inverted Annular
â– D
¿i ,«,1-^..^^^^

23
Table 2.1 continued.
- - multiple
scatter
single
scatter

24
A^ and AQ are as defined for equation (2.1), and
S° is the source strength.
There are several implicit assumptions expressed in
equations (2.6*) and (2.7). They are that a) the
divergence and width of the source beam in the direction
perpendicular to the imaging plane is small, b) the source
is a point body, c) the detector is a point body, and d) the
effects of the detector electronics and MCA are neglected.
Careful construction of the source collimator can insure
that a) is true. Including the effects of b), c), and d),
poses no serious analytic hurdles as far as the forward
problem of predicting the detector output given the density
distribution. However, these details have serious
implications on the resolution obtainable when trying to
reconstruct the density distribution from the detector
output.
,,Equation (2.4) provides a one-to-one correspondence
between the energy of a detected singly scattered photon and
the point of scattering. The relationship is true only if
the source and detector are point bodies. If the so'^ce ana
detector have finite extent, then there is some uncertainty
associated with the scattering angle as depicted in
Figure 2.4. This means that photons scattering at a point
in the density field do not arrive at the detector with a
single energy but with an energy over a small range,
AE = (3E/3©)A0

25
Figure 2.4:
Illustration of the scattering angle blur
problem arising from the use of non-point
source and detector.

26
This "blurring" effect is analogous to taking a picture with
an out-of-focus camera. The amount of uncertainty in the
scattering angle is not constant over the entire density
field nor even along an isogonic line. Consequently, the
blurring cannot be expressed as the convolution of some
spread function with the spectrum and cannot oe removed by
a straignr forward deconvolutional procedure.
There is also blurring of the detected photon spectrum
by the detector and electronics. Unlike the blurring
discussed above, this blurring is generally energy invariant
and can be removed by a deconvolutional procedure. Methods
for removing the detector system blur and correcting for the
energy dependent efficiency are well developed and it is
assumed that the measured data are automatically corrected
for these effects.
Another aspect of the data collection process which
affects image quality is the use of a multi-channel analyzer
(MCA) to record the photon spectrum. The impact of collec¬
ting the spectrum data in an MCA is that the spectrum is not
available as a continuous function but as a discrete set of
values. That is, D(E) becomes the set {Dc;c=l...I} where
Dc = / D(E)dE.
AEc
It is now convenient to introduce the isogonic segment as
shown in Figure 2.5. An isogonic segment is the region
bounded by the isogonic curves which correspond to the

Y
Figure 2.5:
Illustration of an isogonic segment. (Reproduced
by permission. Copyright ® ISA 1985 [21].)

28
energy limits of a single channel in the MCA. Any features
of the spectrum smaller than the width of a single MCA
channel are smoothed over and cannot be recovered.
Generally this does not pose a problem since modern equip¬
ment can collect data with very narrow channel widths.
Parallel Beam Source/Uncollimated Detector
The analysis so far has dealt exclusively with a small
source and a diverging beam. It is relevant to extend the
analysis to parallel source beams since work has been done
using parallel beam illumination with neutron scattering
[20]. Figure 2.6 illustrates one possible way of
constructing a parallel beam source from an array of
narrow beam sources. Alternatively, a parallel beam
source can be seen as the result of moving a wide beam
point source infinitely far away from the pipe. As the
source is moved away, the source-to-detector distance goes
to infinity and by equation (2.5) the radius of an
isogonic curve also goes to infinity. The curvature (i.e.
the reciprocal of the radius) goes to zero and the
isogonic curves become lines or rays which originate at
the detector. The exit paths of the scattered photons (or
neutrons) are, therefore, colinear with the isogonic
lines (see Appendix II). The same equations describing the
detector response for the wide beam source/uncollimated
detector case are valid provided the appropriate limits on
the integrals are used.

29
Discrete Forward Equation
This section deals with the practical aspects of
expressing equation (2.7) in a form which enables computing
the integral of equation (2.6*) using a digital computer.
The details of how the computations are done are presented
in Appendix III. The purpose hereis to show the relation¬
ship between the discrete equations, which are used in
reconstructing an image of a flow cross section, and the
"exact" integral form.
The measured data are available as a discrete set from
which the singly scattered flux spectrum must be extracted.
The extraction process of removing the multiply scattered
flux component and correcting for the detector system
efficiency can be expressed as
sc = where is the number of photons that arrived at the
detector having suffered only a single Compton scattering
collision. The notation denotes a vector with components
S^. The object here is to be able to calculate given the
density distribution in the illuminated cross section of the
pipe. An expression for calculating is derived using
equations (2.7) and (2.7*)
= S °f dEj VnpA^dl .
AEC ME)
(2.9)

30
Beam
Figure 2.6: A parallel beam illumination/uncollimated
detector gamma-ray scattering system.
Isogonic Segment
Figure 2.7: An n x n pixel grid is overlayed the flow cross
section in order to facilitate calculation or
the singly scattered photon flux. Individual
pixels are designated 2^Y a single number k
ranging from 1 to K(=n ).

31
The first step in numerically evaluating equation (2.9)
is to overlay the illuminated flow cross section with an
n x n pixel grid as shown in Figure 2.7. Equation (2.9) is
rewritten as a sum of integrals over the individual pixels
si = S° l / dEf v pkA±A dl (2.10)
kc aec Me)
where pk(x,y) is equal to p(xfy) for points inside pixel k
and is zero elsewhere. The notation kc indicates that the
sum is to be carried out only over those pixels in the
isogonic segment c. The integrals in equation (2.10) can
be expressed equivalently by
= S°i Jdx/dyywpkAiA0PSFc(x,y) . (2.11)
k
The term PSFc(x,y) is the characteristic function of the
isogonic segment c defined as
PSFc (x'y) 0 if
(x,y)
(x,y)
is inside isogonic segment c
is outside isogonic segment c.
Any real imaging system imposes a lower limit on the
spatial resolution below which variations in the density can
not be seen. If the size of the pixels are about the same
as the spatial resolution limit, then it is reasonable (even
necessary) to assume a constant density within each pixel.

32
The remaining integral in equation (2.11) is approximated by
Sc = s° I PkAckPSFck ' (2'12>
k
The factor AQk is an average quantity defined by
Ack
//dxdyuw(x,y)A±(x,y)AQ(x,y)
J/dxdy
k
(2.13)
where the integrals are carried out only over the area of
pixel k. The function PSFck is the fractional area of pixel
k which is inside isogonic segment c or
PSFck = //dxdy PSFc(x,y). (2.14)
k
An implicit assumption in equation (2.12) is that the
attenuation factors do not vary significantly across a
pixel. This assumption is valid except for pixels near an
edge in the density distribution and simultaneously on the
edge of an isogonic segment. The assumption also breaks
down when the size of the pixels are on the order of a
photon mean-free-path. In practice there are no sharp edges
to the isogonic segments due to the scattering angle blur.
The scattering angle blur causes the isogonic segments to
"feather" off at the edges (and overlap) as shown in

33
Figure 2.8. The effect is that the assumption works in real
applications where non-point source and detectors are used.
The definitions for A^ and PSFC^ must be modified to
account for the scattering angle blur. The new definition
for Acj< specifies that the attenuation factors must not only
be averaged over the area of the pixel but also over the
range of possible energies with which the photons can reach
the detector. The new definition of PSFC^ is the fraction
of photons reaching the detector from pixel k with
energies in the range of channel c multiplied by the pixel
area.
Equation (2.12) gives an expression for calculating the
singly scattered photon spectrum for a detector looking at
the scattered photon flux of a known density distribution
illuminated by a wide beam mono-energetic gamma-ray source.
Each channel of the singly scattered photon spectrum is seen
to represent an integral measure of the density distribution
within the corresponding isogonic segment. The goal of the
/
remainder of this report is to do the inverse operation of
extracting information about an unknown density distribution
given measurements of the singly scattered photon flux.
Chapter 3 tackles this inverse problem by correlating the
size and shape of spectra with the different possible flow
regimes. Chapter 4 goes one step further in solving the
inverse problem and looks at ways of inverting
equation (2.12).

34
PSFch(x,y)
y
The scattering angle blur causes the isogonic
segments to feather off and overlap. PSFch(x,y)
is the fraction of those photons which scatter
at point (x,y) toward the detector that have an
energy in the spectrum channel ch.
Figure 2.8:

CHAPTER 3
METHOD OF SPECTRAL MOMENTS
The basic premise of Compton scatter tomography is that
the singly scattered photon spectrum contains information
about the density distribution. If this is true then each
flow pattern should impart some unique characteristics to
the spectrum. Examination of the first few moments of the
spectrum is a simple pattern recognition technique employed
in this chapter to extract these unique characteristics.
Description of the Method
Consider a Compton scattering setup using a wide beam
source and uncollimated detector as illustrated in
Figure 3.1. Mono-energetic gamma rays leave the source and
undergo Compton scattering interactions inside the pipe.
Some of the scattered gammas are received by the detector
and the energy spectrum of the detected gamma rays is accu¬
mulated in an MCA. The total count in each MCA channel is
proportional to the density of the material in the isogonic
segment for the channel. Note that the source beam is
stopped down slightly in order to reduce the amount of the
pipe wall which is illuminated.
The typical two-phase fluid problem consists of a fluid
and its vapor. The vapor is, in general, considerably less
35

36
Figure 3.1: Compton scatter densitometry configuration.
Scattered radiation from the illumination beam
is received by the detector. (Reproduced by
permission. Copyright ® Instrument Society of
America (ISA) 1985 [21].)

37
dense than the liquid. Therefore, an MCA channel represents
an integral measure of the liquid in the isogonic segment.
The singly scattered photon spectrum is a projection of the
liquid distribution onto the energy axis. The energy axis
can, by suitable transformation, be related to any spatial
axis intersecting the isogonic lines.
A convenient reference axis perpendicular to the
isogonic lines is the Y-axis. The transformation relating
the energy axis to the Y-axis is available from
equation (2.5)
Y(E) = r(0(E))+yc(6(E)) . (3.1)
Note that each detector has a characteristic Y-axis.
The fluid distribution in the pipe can be inferred from
the one dimensional image (i.e. projection) by examining the
moments of the spectrum with respect to the Y-axis. The
zeroth moment is a measure of the amount of liquid phase
present. The first moment is a measure of where the center-
of-mass of the liquid phase is. The standard deviation is a
measure of how the liquid phase region is distributed about
the center-of-mass. The moments are not direct indicators
of their respective quantities because of the attenuation
factors.
The zeroth moment is calculated in a straight forward
manner as

38
= IS¿
c
(3.2)
where is a channel of the singly scattered photon
spectrum [see equation (2.8)]. The first moment and
standard deviation are calculated according to the formulae
and

c AY
YdY
c Yc+1
(3.3a)
l si

c AY,
fY,
Y - )2dY
â– c+1
(3.3b)
I si
Figure 3.2 graphically depicts the zeroth moment, first
moment, and standard deviation for annular and inverted
annular flow at 50% void fraction. The two peaks in the
spectrum for the annular flow correspond to the liquid in
the illuminated portions of the annulus separated by the low
density vapor core region. The peak heights are not the
same because the source beam is attenuated before it reaches
the far side of the annulus which results in fewer scat¬
tering collisions taking place there. Geometric attenuation

RELATIVE DETECTOR COUNT
39
Figure 3.2: Spectral moments of the singly scattered photon
spectrum. a) Spectrum for an annular flow;
b) Spectrum for an inverted annular flow.
(Reproduced by permission. Copyright © ISA
1985 [21].)

40
1.00
o
o
cc
o
H
O
Ui
H
UJ
O
Ui
>
H
<
0.75-
0.50
0.25
UJ
CC
< yO > -
0.340 reí. units
r1> \
I'
,0* | \
aA I.
I
I,
i \
i , V
10 15
Y (cm)
20
(b)
Figure 3.2—continued

41
of the photon fluence (i.e. photons/area) for points further
from the source also occurs. Similarly, these same attenua¬
tion terms are applicable with respect to the detector.
Only a single peak for the liquid core region is seen in the
inverted annular spectrum. The peak is not symmetric due to
attenuation as the gamma rays pass through the liquid and to
geometric attenuation.
The zeroth moments are not identical even though the
void fraction is the same. All of the liquid core region is
illuminated in the inverted annular flow. Not all of the
liquid annulus in the annular flow is illuminated. In
addition, the different liquid distribution leads to differ¬
ences in attenuation, solid angle, and differential
scattering cross section. Consequently, the zeroth moment
of the inverted annular flow is larger than that of the
annular flow.
Both the annular and inverted annular flows are
symmetrically distributed within the pipe and so, the first
moment would be expected to lie at the center of the pipe.
However, due to the attenuation effects the first moment is
skewed closer to the source and detector. The standard
deviation of the annular spectrum is larger than that of the
inverted annular. This is expected since the annular
spectrum consists of two widely separated peaks while the
inverted annular spectrum is a more compact single peak.

42
Procedure and Results
Monte Carlo transport calculations are carried out to
test whether flow regime identification is possible by
examining the first few moments of the singly scattered
gamma ray spectrum. Several different flow regimes are
analyzed; in each case the energy spectrum from single
scattering interactions with liquid or vapor inside a pipe
is calculated for two detectors. The detectors are assumed
to be placed symmetrically with respect to the pipe and
source as shown in Figure 3.3. The source-pipe-detector
geometry is the same for all cases.
The source and detector are treated as points for
purposes of calculating the scattering angle. However, for
purposes of calculating the solid angle view of each
interaction, the detector has a radius of 0.5 cm. Detected
scattering interactions within the pipe wall are not
considered since the primary focus is to evaluate the gamma
rays coming from the fluid. However, attenuation due to the
pipe is taken into account.
The flow regimes considered are annular, inverted
annular, stratified, bubbly, and mist. The annular and
inverted annular flows are modeled as a liquid/vapor annulus
and a vapor/liquid central core concentric with the center
of the pipe, respectively. The stratified flow is modeled
as a pipe half-filled with liquid as illustrated in
Figure 3.4. The bubbly and mist flows are simulated by
choosing fifteen randomly sized (0.85
43
bubbles of vapor/liquid in a matrix of liquid/vapor, respec¬
tively. New bubbles are periodically chosen during a cal¬
culation to simulate stochastic flow conditions. The void
fraction is fixed at 50% for the annular, inverted annular,
and stratified flow. Due to the random nature of choosing
the bubbles, the bubbly flow has a void fraction of ~45% and
the mist flow has a void fraction of ~58%. Four calcula¬
tions are made for each flow regime where the width of the
source beam in the plane of the pipe cross section varies
from 4.5 to 28 degrees. The results of the moment
calculations are tabulated in Table 3.1.
Interpreting the Moments
The moments are functions of the magnitude, shape, and
width of the spectra. The spectra change as the source beam
width varies since the amount of the flow cross section
illuminated changes. Stating general trends of the spectra
based on the changing amount of illuminated fluid is diffi¬
cult due to the effects of attenuation. However, the
following two observations can be made. First, the magni¬
tude of the spectrum from a narrow source beam is smaller
than that from a wider source beam. Assuming the source
strength remains constant, the number of photons illumina¬
ting the pipe is less for the narrow source beam. The
result is that the zeroth moment becomes smaller as the
source beam narrows.
Second, the width of the spectrum decreases as the
source width is narrowed. As illustrated in Figure 3.5, the

44
Source
Figure 3.3: One source, two symmetrically positioned
detector configuration used for the Monte Carlo
studies presented in this chapter. (Reproduced
by permission. Copyright ® ISA 1985 [21].)
Figure
Vapor
Liquid
3.4: Source-detector configuration relative to the
stratified flow regime case. (Reproduced by
permission. Copyright ® ISA 1985 [21].)

TABLE 3.1
Zeroth moment, first moment, and standard deviation of the singly scattered photon spectra
for the five flow regimes listed.
Source Beam Width
2 8°
14°
8°
4.5°
Flow Regime

[relative]

[cm]
[cm2 ]
A
0
V





A
o
>i
V

A
Q
M
ANNULAR
Detector 1
0.311
15
12
0.0935
15
15
0.0485
15
14
0.0267
15
14
Detector 2
0.313
15
12
0.0935
15
15
0.0484
15
14
0.0267
15
14
INV. ANNULAR
Detector 1
0.337
16
3.7
0.194
16
3.7
0.108
16
3.5
0.0601
16
3.5
Detector 2
0.340
16
3.7
0.193
16
3.6
0.108
16
3.6
0.0600
16
3.5
STRATIFIED
Detector 1
0.621
17
6.1
0.360
17
6.6
0.207
16
6.9
0.116
16
6.7
Detector 2
1.00
14
5.4
0.353
15
6.3
0.158
16
6.6
0.0789
16
6.4
BUBBLY
Detector 1
0.434
15
8.4
0.204
16
7.6
0.108
16
8.1
0.0597
16
8.2
Detector 2
0.436
15
8.4
0.205
16
7.5
0.107
16
8.0
0.0602
16
8.2
MIST
Detector 1
0.431
16
6.6
0.183
16
6.5
0.108
16
5.9
0.0614
16
5.3
Detector 2
0.423
16
6.9
0.184
16
6.7
0.112
16
5.7
0.0600
16
5.2
Note: The spectra are from a
Cs-137
source
illuminating
í a 12 i
cm diameter
pipe with a 0
.75
cm thick
steel wall
. The beam
divergence
normal
to the
flow
cross
section
is 5
degrees in all cases. The void fraction is nominally 50% for all cases.

46
region within the isogonic segments associated with the
extreme channels of the spectrum is no longer illuminated as
the beam narrows. The effect is more pronounced for the
isogonic segments at low values of the Y-coordinate. If
liquid is present in the affected regions, the result is
that the spectrum narrows and the narrowing can be more
pronounced at the low end of the Y-axis. Such an asymmetric
change in the spectrum width means the first moment shifts
to larger values of Y as the beam narrows. A narrower
spectrum would also be expected to yield a smaller standard
deviation. The magnitude of change in the first moment or
standard deviation, or whether any change occurs at all,
depends on any accompanying change in the shape of the
spectrum. If there is no liquid in the affected regions
then, of course, there is no change in the width of the
spectrum.
A Detailed Example ''
The following is a discussion in physical terms of how
the moments vary with changes in the source beam width. So
that sufficient attention can be given to the changes in
each moment, the discussion is confined to only the strati¬
fied flow regime. Enough insight should be gained to
interpret the data in Table 3.1 for the remaining flow
regimes.
Figures 3.6 and 3.7 are the spectra for the upper and
lower detector, respectively, viewing a stratified flow.
These figures show the progressive changes in the magnitude,

47
shape, and width of the singly scattered photon spectrum as
the source beam width narrows. The magnitudes of the spectra
decrease as the beam narrows because there are fewer photons
illuminating the flow. The magnitudes of the two figures
should not be compared since the vertical scales are diffe¬
rent. Instead the figures should be compared on the basis
of the zeroth moments.
The zeroth moment of the lower detector's spectrum is
larger when the source beam is wide, while the zeroth moment
of the upper detector's spectrum is larger when the source
beam is narrow. The behavior of the zeroth moments can be
explained by considering the effect of attenuation on the
scattered gamma rays and the differences in solid angle each
detector views of the scattering interactions. When the
source beam is widest the average attenuation of the
scattered gamma rays reaching each detector is about the
same. The average position of scattering occurs in the
middle of the liquid and the scattered gamma rays must pass
through almost the same amount of liquid and steel to reach
either detector. When the source beam is narrowed the
average scattering interaction takes place toward the center
of the pipe and the scattered gamma rays must pass through
more liquid to reach the lower detector. Thus, attenuation
of the scattered gamma rays explains why the zeroth moment
of the upper detector is larger when the source beam is
narrow. When the source beam is wide the zeroth moments
should be the same since the average attenuation is the

48
Figure 3.5: Darkened regions are those isogonic segments
which are eclipsed when the more restricted
beam illumination is used. (Reproduced by
permission. Copyright @ ISA 1985 [21].)

RELATIVE DETECTOR COUNT
49
Figure 3.6: Detected photon energy spectra of once-
scattered gamma rays for upper detector viewing
a stratified flow as in Figure 3.4. (Reproduced
by permission. Copyright ® ISA 1985 [21].)

RELATIVE DETECTOR COUNT
50
Figure 3.7: Detected photon energy spectra of once-
scattered gamma rays for the lower detector
viewing a stratified flow as in Figure 3.4.
(Reproduced by permission. Copyright © ISA
1985 [21] . )

51
same. The reason they are not is because the upper detector
is further away and presents a smaller solid angle to each
scattering interaction.
The behavior of the first moments can be explained in
the following manner. Observe that the spectrum for the
widest source beam is peaked at about 13 cm for the lower
detector and at about 17 cm for the upper detector. Due to
the source-pipe-detector geometry the isogonic segments do
not all intersect the liquid region. The projection of the
liquid region onto the corresponding Y-axis for the two
detectors is not symmetric. As the source beam narrows both
spectra become more uniform since approximately the same
amount of liquid is illuminated in each isogonic segment.
The remaining slight peak is due to the non-uniform amount
of attenuation and to solid angle considerations. The more
uniform a spectrum the more the center-of-mass coincides
with the center of the spectrum. The result is that the
first moment shifts to higher Y values for the lower
detector and shifts to lower Y values for the upper
detector. Note also the preferential narrowing of the
spectra which tends to enhance the changes in the first
moments.
The standard deviation for both detectors' spectra
increases uniformly as the source beam narrows from 28 to 8
degrees and then decrease slightly as the beam narrows
further to 4.5 degrees. When the source beam is widest both
spectra are peaked near the center-of-mass which gives a

52
relatively small standard deviation. The standard deviation
increases as the source beam narrows because the spectra
become more uniform. The slight narrowing of the spectra
leads to a smaller standard deviation. The two competing
effects result in the observed behavior.
Using the Moments
The purpose of calculating the moments is to identify
the flow regime automatically and without any a priori
knowledge. Inter-comparing the moments of all the flow
regimes at each particular source beam width suggests that
using the widest source beam offers the greatest discrimina¬
tion between flow regimes. A close agreement between the
moments of the upper and lower detectors is seen in all the
data for the annular, inverted annular, bubbly and mist
flows where the fluid distribution is symmetric. The
stratified flow is asymmetric, but only for the wider source
beams (i.e. 14 degrees) is the asymmetry discernible from
comparing the moments of the upper and lower detector. The
annular flow is discernible from the other symmetric flows
with the wide source illumination because of its large
standard deviation. The inverted annular flow conversely is
distinguishable because of its small standard deviation.
Discriminating between the bubbly and mist flows is not
possible on the basis of the first moment or standard
deviation due to their similar nature. Bubbly flow is
expected when the void fraction is low (i.e. < 50%) while
mist is expected when the void fraction is high. The zeroth

53
moment is sensitive to the void fraction and could be used
to determine which flow exists.
The results presented so far show that flow regime
identification can be made using the singly scattered photon
spectra from two detectors sensing the emerging scattered
photon flux from a flow cross section illuminated by a wide
beam source. A tomographic image of the flow can be
produced using a simple model of the flow cross section for
each flow regime. The consistency between the actual flow
regime and the flow regime as identified by the spectral
moments can be checked using the solution to the forward
problem [i.e. equation (2.12)]. The calculated singly
scattered photon spectra for the two detectors based on the
identified flow regime should match the measured spectra.
If the calculated and measured spectra are not close enough,
then the flow cross section model is adjusted accordingly.
For example, if the spectral moments indicate that the flow
is annular, then the size of the annulus is adjusted until
the measured and calculated spectra match. The next chapter
develops other methods of producing such tomographic images
of the flow. These other methods are more systematic in
adjusting the flow cross section image so that the measured
and calculated singly scattered photon spectra match.

CHAPTER 4
THE RECONSTRUCTION PROBLEM
A general statement of the problem at hand is
Given a set of singly scattered photon flux measure¬
ments, {S1}, made using a wide-beam source and
uncollimated detectors, reconstruct the liquid-vapor
phase distribution {p} that produced {S^.
Ideally the set of flux measurements is acquired with
a single exposure. That is, the source and detector
locations are fixed and there are no moving parts. Insight
into this problem can be gained by considering the following
formulation of the forward problem from equations (2.6*)
and (2.7):
Sl=[T]p
(4.1)
where [T] is the matrix which mathematically describes the
scattering transport process, §■*■ is the vector of singly
scattered photon flux measurements, and p is the vector of
pixel densities. The obvious solution of inverting [T],
that is,
p=[T] 1 S1
(4.2)
54

55
is not strictly possible since the elements of [T] are
functions of p. The problem at hand is nonlinear. Until
now the general approach has been to avoid dealing directly
with the nonlinearity of [T],
One often-proposed method (e.g. [24, 25]) is to use
multiple source-detector locations and/or source energies in
order to cancel out the nonlinear aspects of [T]. There are
numerous practical limitations in implementing this approach
for dynamic measurements. Another proposed method is to
assume a density distribution based on prior experience and
use this in computing [T]. Since [T] is a function of the
density, a bias in p is likely for the high contrast
reconstructions expected from two-phase flows. A charac¬
teristic of these methods is the idea of calculating the
density distribution in an algebraic, once-through fashion.
It is proposed here that the problem be viewed as a
tomographic reconstruction problem where the measured data,
s} represent projections of the phase distribution, p, onto
non-physical energy axes. The projection data are
incomplete since only a few detectors are employed. Simila¬
rly the elements of [T] are at best known only approximately
since they depend on the unknown phase distribution p. When
the available information is incomplete and/or imperfect, an
algebraic solution may not exist or several solutions may
exist. These are familiar situations in CT and the analytic
techniques are already developed to handle them. The analo¬
gous projection operator, [T], of current CT applications is

56
mathematically well-behaved (i.e. linear, single-valued,
etc.). The CT reconstruction methods must, therefore, be
modified to handle the nonlinearities of Compton scatter
tomography (CST). Several CT reconstruction techniques are
adapted to CST in this chapter.
General Solution Methods
The basic methods of CT for solving the reconstruction
problem can be grouped into five categories:
1. Fourier transform: exponential Radon transform [26],
filtered backprojection [27].
2. Direct scanning [12, 14, 15].
3. Matrix inversion: successive approximation [20],
maximum entropy (see Appendix IV) [28], generalized
inverse, and pseudo inverse [29].
4. Series expansion [30]: algebraic reconstruction
technique (ART), simultaneous iterative reconstruc¬
tion technique (SIRT), and iterative least squares
(ILS).
5. Statistical: Monte-Carlo backprojection.
The transform techniques are recognized as being the
more elegant, mathematical reconstruction methods. Trans¬
form methods in general require a large number of projec¬
tions and do not show good noise immunity. The cost of
detectors and the physical limitations in mounting and
cooling the detectors imply a limit on the number of projec¬
tions available. Practical limits on source size, limits on
detector size, and desired framing rate are all going to

57
result in statistical noise in the measurements. These
features of using Compton scatter to image dynamic two-phase
systems make transform methods unsuitable as a recon¬
struction technique.
Direct scanning techniques inherently require a narrow
collimated source beam and a moving source detector arrange¬
ment. A narrow collimated source is less photon efficient
than a wide collimated source. This inefficiency coupled
with the time necessary to translate the source and detec¬
tors means that a direct scanning system is slower than a
comparable wide-beam source, uncollimated detector system.
Since one goal is to measure dynamic two-phase flow, the
faster measurement system is the obvious choice. There is
another reason for dismissing direct scanning techniques
that involves accounting for the multiply scattered photon
flux. The structure of the flow regime is more strongly
imprinted on the multiply scattered photon flux when a
narrow-beam source is used than when a wide-beam source is
used. Any structure in the multiply scattered photon
spectrum means that approximating the multiply scattered
photon spectrum by a general function B(E) [see equation
(2.6*)] is not possible. Removing the multiple scatter
component in direct scanning systems, therefore, requires
more computation and is more prone to error.
Matrix inversion and series expansion reconstruction
techniques are equivalent in that a solution is sought by
iteratively improving an initial guess of p (and hence [T])

58
until some integral transformation of the reconstructed
phase distribution and the true phase distribution are
judged close enough with respect to the measurement error.
The difference between matrix methods and series expansion
methods is that at each iterative step matrix methods update
all the pixel values at once, whereas series expansion
methods update the pixel values sequentially. Matrix
methods usually exhibit better convergence rates than series
expansion methods since they can take into account higher
order effects in a more straightforward manner. Techniques
from either of these categories appear to be well suited for
use in a CST system.
A Closer Look
This section details how specific CT techniques are
adapted for use in CST. The first technique to be
considered is backprojection. Generally the initial step in
a reconstruction algorithm is a backprojection operation
which yields a very smooth representation of the true image.
Much of the fundamental differences between traditional CT
and CST are evident in the way the projection data are
handled in order to perform a basic backprojection opera¬
tion. The second technique considered is Monte-Carlo
backprojection (an original concept by this author and Dr.
E. E. Carroll at University of Florida). The result of
Monte-Carlo backprojection is a higher contrast image than
obtained with traditional backprojection.

59
The third part of this section is a general discussion
of matrix inversion techniques. The purpose of this
discussion is to introduce the basic iterative nature of any
practical solution strategy imposed by the nonlinearity of
CST. The fourth part of this section adapts three series
expansion techniques (ART, SIRT, and ILS) for use in CST.
Backprojection
Consider a projection formed by the conventional CT
scanner shown in Figure 4.1. The pencil beam source irra¬
diates the object and the transmitted radiation is measured
at points in the image plane. The source and detector are
moved in the direction indicated to scan the entire object
and form a single projection. Other projections are formed
by rotating the source and detector around the object and
taking another scan. The number of photons reaching the
detector at each point in a projection is a function of the
average optical path through the object. The functional
relationship is
Ijc = IQexp{-/ p(s)vit(s)ds} (4.3)
£jc
where I^c is the output of the detector at position c in
projection j,
IQ is the number of photons leaving the source,
p(s) is the density at point s,
Pt(s) is the total mass attenuation coefficient of the
material at s, and
indicates that the integral is to be carried out

60
Figure 4.1: Example of a conventional transmission CT
scanner. The source and detector move parallel
in the direction indicated by the straight
arrows, and also rotate about the object being
imaged.

61
along the line between the source and detector
when at position c in projection j.
In practice one does not work with I_:_, but instead with R â– 
J
given by
K
Rjc = -ln(Ijc/I0) = yt l PkAsjck (4.4)
k-1
where pk is the density in pixel k,
As^ck is the length of line il^c contained in pixel k,
and
K is the total number of pixels.
The value of ASjck is nonzero only if line &jc crosses pixel
k and so the sum in equation (4.4) is equivalent to the
integral in equation (4.3). Since the liquid and vapor are
of the same material, the mass attenuation coefficient is a
constant and is taken outside the integral. The average
density, p^„, along line SLj.„ can be computed from equation
(4.4) as
Jc
= R
DC
K
/{J-i-f- I Asjck}
k=l
^jc^kjc
(4.5)
A backprojection image is made by forming a composite of the
average line densities, p.;„,
J ^

62
J C ( j )
bPk =i l l
j=l c=l
PjcAsjck
}/
J
I
C( j )
l As_i
j=l c=l
jck
(4.6)
where bpk is the density of pixel k in the reconstruction,
and ASjc]< is used as a weighting factor.
An analogous relationship for computing a backpro-
jection image in CST can also be developed. The quantity
corresponding to R^_ in CST is the singly scattered photon
flux values, S^, given by
K
s°j MPkAjkPSFjckAs
JC — -L
(4.7)
where u
S°
Pk
As
Ajk
PSF^ck is uieuractiun ui me pnuuona arrxvxng at
detector j from pixel k which have an energy in
the range of channel c (Ec to
Note that PSFjck is zero if it is not physically possible
for a singly scattered photon from pixel k to arrive at
detector j with an energy in the range Ec to Ec+p So, the
sum of equation (4.8) represents an integral over the
isogonic segment c of detector or projection j and is
equivalent to the sum of equation (4.5). Another point to
is the incoherent scattering coefficient,
is the uninteracted source flux,
is the density in pixel k,
is the length of one side of a (square) pixel,
is an average quantity which takes into account
all aspects of photon transport from the source to
pixel k, and from pixel k to detector j, and
+» U -P v" o /■> 4» i y> y* -P 4» K tr\ 4» j*» v\ e~*

63
note about equation (4.7) is the that double subscript jc
can be replaced by an equivalent single subscript g,
sg = s°J PkAgkPSFgkis (4'7
}C — J.
which is convenient in formulating the problem in matrix
notation as in equation (4.1).
The average density in an isogonic segment, Pjc, is
computed as
pjc
=
K
Dc
/ASS° 1 iAjkPSFjck
(4.8)
and similarly the backprojection image densities, bp^, for
CST are computed as
bPk =
J C ( j )
I I P
1=1 c=l
jcPSFjck
J c { j )
l l PSF
j=l c=l
jck
(4.9)
A basic problem with equation (4.8) is that the density
distribution must be known in order to calculate the Ajj.
factors. A reasonable solution is to use a homogeneous
distribution with an average void fraction that is either
assumed or gained from an independent measurement. The
average void fraction of a CST backprojection image must

64
then be checked for consistency with the initial void frac¬
tion. Any inconsistency can be resolved by some iterative
refinement. Equations (4.5) and (4.6), in contrast, require
no knowledge of the average void fraction and a one-step
backprojection computation is possible in CT.
Monte-Carlo Backprojection
The motivation for using a statistical approach to form
a backprojection image comes from three observations.
First, the singly scattered photon flux, Sjc, is directly
proportional to the average density in the isogonic segment
jc. Second, any detected photons are due almost entirely to
scattering in the liquid phase since the density ratio of
liquid-to-vapor phase for water is usually 1000 to 1 in
nuclear reactors. Third, the features of one detected
spectrum can be cross correlated to the features of another
spectrum detected at another location. These observations
are summarized by the following example. If a single peak
is observed in channel c of detector j and a similar peak is
observed in channel c' of detector j', then the liquid phase
is localized to the region where isogonic segment jc
intersects segment j'c'.
Monte-Carlo backprojection treats each detected
spectrum as an independent one dimensional probability
density function of where the liquid phase is located. A
reconstruction of the liquid phase distribution is made by
statistically cross correlating the detector spectra. The
difference between Monte-Carlo backprojection and

65
traditional backprojection is that traditional backprojec-
tion ignores the cross correlation information.
The initial step in Monte-Carlo backprojection is to
form an independent probability density function (pdf) from
the modified average isogonic segment densities of each
detector
pdf j (c)
Dc
C ( j )
I P
c=l
DC
(4.10)
The average isogonic segment density is modified in the
following manner:
, K
S±/ AsS° l A4
J k=l
jk
(4.8*)
in order to account for the non-uniform spatial width of the
isogonic segments.
Assume that a fine-mesh square pixel grid is employed
and that a pixel is either liquid filled or vapor filled.
Provided a fine enough mesh is used, most density distribu¬
tions can be reasonably represented by such a scheme. The
number of pixels, N, which are liquid filled is easily
computed from an estimate of the average void fraction, a,
N = round((l - a)*K)
(4.11)

66
where K is the total number of pixels. All the pixels are
assumed vapor filled to start and the determination of which
N pixels to designate as fluid filled proceeds as follows:
1.Given a set of J random numbers {r } , an associated
set of channel numbers {ch-} is generated by
inverting the pdfj's
chi
r^ á £J pdf.; (c) ; j = l, 2, J (4.12)
J c=l J
There are J detectors and one channel from each
detected spectrum is chosen by equation (4.12).
2.Associated with each channel chosen is a set of
pixels, Set(ch'), which contains all the pixels
within the corresponding isogonic segment. A new
set is formed by taking the intersection of all J of
these Set(ch-:)'s. The intersection set, I-Set,
contains a total of M pixels. Taking a new random
number, rJ+^, the mth element
m = integer(M*rJ+1) + 1 (4.14)
of I-Set is chosen as the pixel to be fluid filled.
If M is zero, then no pixel is designated as being
fluid filled.
3.Steps 1 and 2 are repeated sequentially until the
necessary N pixels are turned on.
Figure 4.2 compares the results of Monte-Carlo and
traditional backprojection. The Monte-Carlo reconstruction
is not as blurred and better represents the actual flow
regime. Although the reconstructed image resembles the
correct flow regime, a comparison of the calculated and
detectedt singly scattered photon flux spectra shows that
there is still room for improvement (Figure 4.3).

67
(b)
(c)
Figure 4.
^.Detector
--Detector
2: Reconstruction of a 50% void
flow regime. a) Monte-Carlo
b) Simple backprojection; c)
locations.
fraction annular
backprojection;
source and detector

68
Figure 4.3: Comparison of the calculated singly scattered
photon flux for the Monte-Carlo backprojection
reconstruction of Figure 4.2 to the "measured"
spectra, a) Upper detector; b) Lower detector.

69
The addition of a feedback mechanism to the pixel
selection process is one possible improvement. Instead of
randomly choosing a pixel from the intersection set, those
pixels which improve the calculated spectra could be chosen
preferentially. A more complete evaluation of the Monte-
Carlo reconstruction method and possible variations is left
for future work.
Matrix Methods
Generally matrix inversion methods are avoided in CT
due to the large size (i.e. > 1 000 x 1000 ) of the matrices
involved. Even in CST where there are only a few
projections and, with a coarser pixel mesh, the matrices can
be on the order of 200 x 200. Inverting matrices of this
size requires large amounts of storage, long computing
times, and is prone to round off error problems. Whatever
the difficulties in implementing a matrix method solution,
the compact notation is useful.
There are many ways of formulating the reconstruction
problem as a series of matrix operations. The most basic
approach is to treat the problem as an optimization problem.
A search is made for the density distribution, p*, which
minimizes, either directly or indirectly, the difference
between the measured data, S^, and the calculated singly
scattered photon flux, F, that is,
All data in this report, unless otherwise stated, are
the result of Monte-Carlo photon transport calculations
(i.e. numerical experiments) and not physical
measurements.
t

70
min |3^ - F|
P
or min |- [T(p)]p|
P
(4.15a)
subject to the constraints
Pvapor < Pk < Pfluid
k=l, 2
/ • • • t
K.
(4.15b)
The notation |..| represents some norm functional.
Any reconstruction algorithm in CST has a structure
similar to that depicted in Figure 4.4. Note that the
nonlinearity of the problem requires the recalculation of
the transport matrix for each new density distribution, pm.
Changes in the density distribution directly affect the
attenuation of the source photons and scattered photons
which is accounted for in the transport matrix. If the
variation in pm is small so that [T] is approximately
constant, then it may be possible to skip the recalculation
step for a few iterations. In contrast, the corresponding
matrix in CT is only a function of the system geometry and
remains constant throughout the calculation.
Another point at which the nonlinearity of CST needs to
be addressed is in the calculation of Spm+^. One method of
calculating Aprn+'*- is to use the method of steepest descent
5pm+1 = (Jacobian of F(pm))~1•(S1 - F(pm)) . (4.16)
The transport matrix can be treated as a constant or as a
function of the density in the calculation of the Jacobian.

71
Figure 4.4: Flow chart outlining the general structure of a
reconstruction algorithm for Compton scatter
tomography. The superscript m denotes the
iteration index.

72
Treating [T] as independent of p means a smaller step size,
f, is necessary and that Ap does not point directly toward
the path of steepest descent. Including the dependence
of [T] on p results in a significant increase in the storage
requirements and in computing time.
Algebraic Reconstruction Technique
ART is historically the first of the series expansion
techniques applied to commercial CT devices [31]. There are
two basic forms of ART algorithms: multiplicative and
addative. Presented here is a variant of the additive
algorithm.
A single ART iterative step involves correcting the
density vector by adding to it a correction vector
-m+! = -m
cm+l•wm+1
(4.17)
where cm+^ is a constant correction factor and wm+-*- is a
weighting vector. The correction factor is computed so that
one of the calculated singly scattered photon flux values is
consistent with the corresponding measured value, that is,
S
1
g
T?nt+1
Fg
I- yrrm+lrim+l
1~ ¿Tgk Pk
K
(4.18)
Subsequent iterations correct the density vector with
respect to the other measured values, one at a time. The
subscript g is incremented cyclically

73
gm+-*- = (gm + 1) modula G + 1
with each iteration in order to cover all G data points.
The correction factor can be derived by substituting
equation (4.17) into (4.18)
It
k
m+lTTm+l
gk wk
The transport matrix elements are unknown until the new
density vector is known which in turn requires the
correction factor. One solution to this predicament is to
use the transport matrix calculated during the previous
iteration, that is,
,m+l _
el _ Yrpin _.m
sg ¿Tgk Pk
I
mm wm+l
Tgkwk
.m+1
(4.19*)
This approach is valid whenever the transport matrix does
not change much from one iteration to the next, as when
close to the final solution.
There are several possible choices of weighting
vectors. The simplest choice is an identity vector,
however, this results in all the pixel values being altered
including pixels which do not affect the sum in equation
(4.18). An identity vector modified such that

74
g = g
m+1
(4.20a)
sensibly corrects only those pixels contained in the
isogonic segment corresponding to S¿. Another reasonable
y
weighting vector is
g = g
m+1
(4.20b)
which effectively weights a pixel according to the proba¬
bility that a scattered photon from that pixel is counted in
Sg. A third possible weighting vector is
m+1
(4.20c)
g=g
which is analogous to the fractional ray length weighting
scheme of CT [32]. The weighting scheme chosen here and
elsewhere in this report is that of equation (4.20b).
Figure 4.5 graphically shows how an ART solution
proceeds for a two data point, two pixel situation. The
lines L^1 and Lâ„¢ correspond to the two simultaneous equations
given by
(4.21)
The lines change as the density vector is corrected since
the transport matrix is a function of the current density
distribution iterate. The lines tend to stabilize as the

1.0
p2
Figure 4.
0.2
U1
J
0.2
_J
0.4
0.6 0.8 1.0
Pi
Graphic demonstration of how an ART solution proceeds for a two data point,
{L^}, two pixel, {pi}, solution. Superscripts identify the iteration number.

76
final solution is approached, indicating that the transport
matrix is only varying slightly.
Simultaneous Iterative Reconstruction Technique
SIRT is similar to ART except that the density vector
is corrected using all the projection data, S’*-, simulta¬
neously. SIRT attempts to combine a complete cycle of ART
iterations (i.e. G iterations) in calculating the correction
vector. The correction vector elements are given by
m+l _
Pk
i
9
tVal,((sg ' f5,/(
i <’V°§>
(4.22)
where at is the measurement variance. The correction vector
y
resembles a backprojection image of the error vector
em = s1 - Fm
(4.23)
which can be seen by comparing equations (4.22) and (4.9).
Reference 30 indicates that some method of damping may
be necessary to avoid a solution which oscillates and fails
to converge. Once the correction vector is obtained it is
multiplied by a constant damping factor, f. The damping
factor is calculated so that the variance weighted euclidean
norm of the error vector is minimized
min | S1 - [Tm](pm + fm+1Apm+1)|. (4.24)
f

77
The solution to equation (4.24) is
(4.25)
g k
The proposed SIRT algorithm for CST is
where the correction vector is given by equation (4.22) and
the damping factor by equation (4.25).
Iterative Least-Squares Technique
The ILS technique computes a correction vector which
minimizes the variance weighted euclidean norm of the error
vector. This can be expressed as
(4.27)
The correction vector which minimizes U is found by setting
the derivative of U with respect to the correction vector
equal to zero. This leads to the set of equations
k'=l
K (4.28)
Solving the system of equations of (4.28) would violate the
principal motivation of going to an iterative approach,

78
which is introduced to avoid the computational overhead and
pitfalls of matrix methods. The usual course of action is
to assume that all the the correction vector elements are
zero except the element corresponding to the current pixel
k\ Equation (4.28) becomes
k '=1
(4.28*)
K
Equation (4.28*) overestimates the magnitude of the indivi¬
dual elements of the correction vector which necessitates
the introduction of a damping factor as is done for SIRT.
The form of the ILS damping factor, which is calculated
according to equation (4.24) also, is
IÍ(£g/ag) (4.29)
9
iKÍT^Ap^h/Ogí2
The proposed ILS algorithm for CST is
Pk+1 = raid ÍPvapor' Pk + fm+1Apg+1, Pfiuid^ (4'30)
where the correction vector is given by equation (4.28*) and
the damping factor by equation (4.29).

79
Comparison of Reconstruction Techniques
The ART, SIRT, and ILS reconstruction techniques are
compared here in order to select one of them to use in a
detailed evaluation of CST. Each of the techniques are used
to reconstruct the two flow regimes shown in Figure 4.6.
The flow regimes are analogous to the "hot spot" and "cold
spot" phantoms used in CT studies. The criterion used to
judge the best techniqué is the euclidean distance between
the true density vector and the reconstructed density
vector. The rate of convergence is considered as a
secondary criterion.
The projection data (spectra) used for each technique
are generated by Monte-Carlo transport calculations. The
singly scattered photon spectra calculated for the two
detectors are shown in Figures 4.7 and 4.8. Table 4.1 lists
the pertinent "experimental" conditions.
The reconstructions of the hot-spot and cold-spot flow
regimes are shown in Figures 4.9 and 4.10, respectively.
All the reconstructions exhibit undesirable streaking
artifacts. The artifacts are due to the lack of edge infor¬
mation from using only two views. The reconstructions of
the hot-spot are better than those of the cold spot. Intui¬
tively this is expected by observing that the signal-to-
noise ratio is higher for the hot-spot spectra. The ART
results have more pronounced streaking artifacts than either
the ILS or SIRT results. The poorer performance of ART is
in keeping with similar CT phantom study tests done for

(dimensions in centimeters)
x
00
o
Figure 4.6: Source and detector configuration used for the hot-spot and cold spot tests.

Energy [keV]
(a)
(b)
oo
The measured singly scattered photon spectra for the hot-spot test,
a) Spectrum from the upper detector; b) Spectrum from the lower detector.
Figure 4.7:

(a) (b)
Figure 4.8: The measured singly scattered photon spectra for the cold-spot test. The
ordinate scale is the same as for Figure 4.7. a) Spectrum from the upper
detector; b) Spectrum from the lower detector.

83
TABLE 4.1
Data pertaining to the hot-spot and cold-spot experiment.
Detectors:
Number of detectors
Detector type
Diameter
Active thickness
MCA channel width
2
planar
1 cm
0.1 cm
1 keV
Source:
Isotope
Cs-137
Physical shape
cylinder
Diameter
0.31 cm
Height
0.63 cm
Material Density: Hot-spot
bubble density
1.155
g/cc
bulk density
0.0
g/cc
Cold-spot
bubble density
0.0
g/cc
bulk density
1.155
g/cc

SIRT
Figure 4
V:«'f
ART
r'*'-
â– s~. â– 
â– f
.9: Reconstructions of the hot-spot flow pattern
by the ILS, SIRT, and ART algorithms.

85
Figure 4.10:
Reconstructions of the cold-spot flow pattern
by the ILS, SIRT, and ART algorithms.

86
medical applications. There is no discernible difference
between the ILS or SIRT results.
The rate of convergence and accuracy of the reconstruc¬
tion methods is shown in Figure 4.11. The discrepancy
function is defined as the euclidean distance between the
true flow regime and the reconstructed flow regime. The
spectra difference function is defined as the euclidean
distance between the measured spectra and the spectra calcu¬
lated using the reconstructed flow regime. Both the discre¬
pancy and spectra difference functions are normalized by
their values for the initial density distribution, which is
an all vapor condition. The convergence of all the recon¬
struction algorithms for the cold spot test is similar. The
ILS and SIRT methods show better convergence properties than
the ART method for the hot spot test as judged by the
discrepancy function. The ILS and SIRT method are more
desirable than the ART method due to their consistent
performance. There is no significant difference in the
performance of the ILS or SIRT methods based on examining
the discrepancy or spectra difference functions. The choice
of whether to use either the ILS method or the SIRT method
appears to be one of personal preference.
This chapter successfully demonstrates that CT recon¬
struction algorithms can be used with CST. As vital as the
reconstruction algorithms are, they are not sufficient for
implementing a CST system. The next chapter considers the
effect of noise in the spectra on the reconstruction process

Discrepancy
o ILS
â–¡ SIRT
A ART
Comparison of the accuracy and convergence of the ILS, SIRT, and ART
algorithms. a) Hot-spot flow regime; b) Cold-spot flow regime.
oo
Figure 4.11:

Figure 4.11—continued.

89
and looks at some of the other practical aspects of
designing a complete CST system. The ILS algorithm is
chosen as the reconstruction algorithm in the discussions of
the next chapter and the remainder of the dissertation.

CHAPTER 5
COMPTON SCATTER TOMOGRAPHY
Three major issues of interest in designing a CST
system or in evaluating the suitability of CST for a
particular purpose are accuracy, image resolution, and image
contrast. The nonlinearity of CST makes it difficult to
quantify any of these for all values of void fraction, flow
regime, position in space, etc. A complete analysis
requires varying all pertinent parameters combinatorially.
The following discussions do not attempt a complete quanti¬
tative analysis, but seek to identify general trends and
limitations of CST.
An Example
In order to illustrate how the non-linear aspects of
CST affect the flow image reconstruction, consider the case
of a small error in the source strength specification. The
detector response is directly proportional to the source
strength [see equation (2.6*)]. If the correct flow distri¬
bution along with too low a value for the source strength is
used in calculating the detector response [see equation
(2.12)] then the calculated spectrum will be uniformly lower
than the measured spectrum (Figure 5.1). The reconstruction
algorithm is only able to compensate for the error in S° by
90

Counts
Correct
Figure 5.1:
Result of specifying too low a source strength to the reconstruction program.
The calculated spectrum is based on the correct flow distribution, but with
S° 10% too low.

92
adjusting the fluid density in the individual pixels.
Examining equations (2.1) and (2.6*) it is seen that the
density enters the detector response equation as a linear
factor and as an argument of the exponential attenuation
factors. If the attenuation factors did not exist then the
low value of S° would result in a uniform increase in the
density throughout the flow image thereby uniformly
increasing the calculated spectrum to match the measured
spectrum. However, the attenuation factors are present and
as the density increases the attenuation increases (i.e. the
magnitudes of the attenuation factors decrease) requiring
more of an increase in the density than if the attenuation
factors are not present. Since attenuation is not uniform
over the flow cross section this leads to a non-uniform
change in the density or, equivalently, a change in the
density distribution. The constraint conditions placed on
the minimum and maximum allowable density values [see equa¬
tion (4.30)] can exacerbate the error in the density distri¬
bution if some of the pixels are constrained from changing
while others are not. Specifying a source strength value in
the input to the reconstruction algorithm which is lower
than actually used in making the measurements results in the
reconstructed density values being too high and the density
distribution being incorrect (Figure 5.2).
Accuracy
The topic of accuracy is confined here to examining the
extent to which errors in the input data result in errors in

93
(a)
(b)
y
x
Figure 5.2
Reconstruction of an annular flow using an incor¬
rect source strength value. a) S° 10% too high;
b) S° correct; c) S° 10% too low.

94
the reconstructed density distribution. Such an approach
views error propagation as a concern in the operation of an
imaging system rather than in its design. The design
process is limited to minimizing the errors in the measured
data due to the lack of analytic tools to handle nonlinear¬
ities. In recognizing this limitation the approach taken
here is to look at a specific system and its response to the
expected types of error.
The specific system chosen is designed to measure the
two-phase distribution of water in a pipe. The system
consists of two detectors and a source placed as shown in
Figure 5.3. No actual experiments are done using this
setup, instead the detector responses are generated by a 3-D
Monte-Carlo photon transport computer code. The various
components of the setup are modeled in the transport calcu¬
lation in the following manner. The source is cesium-137
which has a single gamma ray emission energy of 662 KeV.
The shape of the source is a cylinder 0.64 cm long and 0.32
cm in diameter. A perfectly absorbing shield collimates the
source into a wedge beam which illuminates a thin cross
section of the pipe. The two detectors are 1 cm in diameter
and 0.1 cm thick. These high energy resolution detectors
sense the scattered photon flux. There are a total of 132
channels each 1 keV wide in each spectrum and due to
scattering angle blur it is estimated that this represents
126 independent channels.

y
Figure 5.3: The source-pipe-detector configuration used in evaluating the
effects of flow regime, void fraction, and noise on CST.

96
The presence of the pipe wall is neglected in the
Monte-Carlo calculations and in reconstructing the flow
cross sections. In a real experiment the pipe would be a
fixed scattering (and absorbing) material whose dimensions,
location, and density are known and can be accounted for in
the reconstruction process. Ignoring the pipe reduces the
computer runtime and programming effort considerably.
The most obvious source of error is the statistical
noise in the spectra due to using radiation in the imaging
process. Systematic errors may also be present in the
spectra due to poor instrument calibration, dynamic bias
[33], or spectra processing errors. The only other data
required by the reconstruction algorithm are the physical
dimensions of the system. Since there are no moving parts
in the system, these dimensions can be measured quite accu¬
rately and do not represent a source of significant error.
The system as defined leaves only three operational
parameters which can vary-- data acquisition time, flow
regime, and void fraction. Flow regime and void fraction
are the two quantities being measured while data acquisition
time is the only parameter under control of the operator.
The error in a spectrum is closely related to the data
acquisition time. The shorter the data acquisition time the
fewer the photons detected and the higher the statistical
noise in the spectrum. The longer the data acquisition time
the more the stochastic flow variations introduce dynamic
bias and blur the spectra. The data acquisition time is

97
chosen to balance the statistical and systematic errors in
order to minimize the resulting error in the reconstructed
density distribution.
In the following analysis the error in the recon¬
structed density distribution due to errors in the spectra
under different flow regime and void fraction conditions is
considered. Table 5.1 lists the 3 parameters being varied
and their range of values. The four flow regime models are
illustrated in Figure 5.4 as they should appear when
correctly reconstructed. These four regimes cover a broad
range of two-phase flows and most regimes can be viewed as
combinations or extreme cases of these basic types. Two
particular cases of the stratified regime are considered;
stratification symmetric and stratification asymmetric with
respect to the two detectors. These two cases are chosen in
order to judge the sensitivity of CST to asymmetries in the
flow. Three separate void fraction cases (i.e. 0.3, 0.5,
and 0.7) are looked at for each of the five flow regimes.
Figure 5.5 is an example of the detected scattered
photon spectra for the case of annular flow and void
fraction of 0.5. Similar plots for all fifteen flow regime-
void fraction cases are given in Appendix V. Figure 5.6a is
the reconstructed image of the flow cross section using the
spectra in Figure 5.5. A 15 x 15 pixel array is used and
the number of dots is linearly proportional to the density.
It is known a priori that 6 pixels in each corner are voided
which leaves 201 unknown pixel densities to calculate. The

98
TABLE 5.1
The three parameters
propagation analysis,
being examined as part of the error
Parameter
Range of values
Flow regime
Annular
Inverted annular
Bubbly
Stratified
Void fraction
0.3
0.5
0.7
Spectrum noise
No added noise
10% additive random noise and
25% magnitude decrease
30% additive random noise

99
Figure 5.4: The five flow regime cases at 50% void fraction
as they should appear when correctly recon¬
structed. a) Annular flow; b) Inverted annular
flow; c) Bubbly flow; d) Symmetric stratified
flow; d) Asymmetric stratified flow.

100
Energy [keV]
Figure 5.5: Example of the singly scattered photon spectra
used as input to the reconstruction program.
This particular spectrum is from a 50% void
fraction annular flow.

101
Figure 5.6: Reconstruction of the annular flow regime using
the spectrum in Figure 5.5 with the three noise
conditions. a) No added noise; b) 10% additive
noise and 25% magnitude decrease; c) 30% addi¬
tive noise and no magnitude decrease.

102
ratio of unknowns (pixels) to knowns (spectrum channels) is
approximately 0.80.
The ILS algorithm is employed in all the reconstruc¬
tions in Figure 5.6 and this chapter. The reconstruction
algorithm is an integral part of a CST imaging system and
the sensitivity of CST to the experimental parameters is
somewhat dependent on the algorithm chosen. In the most
general system in which no assumptions are made about the
flow regime the ILS algorithm should perform as well or
better than any other algorithm.
Three different noise or error cases are considered.
The first case is the detected spectra with only the
inherent noise from the Monte-Carlo calculations (less than
6%). The second case adds a normally distributed random
noise component to each measured spectrum and introduces a
systematic error by decreasing the magnitude by 25% (i.e.
source magnitude error). The standard deviation of the
noise component is 10% of the measured spectrum's magnitude.
Similarly, the third noise case adds a 30% normal random
noise component, but introduces no systematic error. The
normally distributed random noise is generated by the method
of adding five uniform random deviates, correcting the
values to have zero mean, and adjusting the standard devia¬
tion to the appropriate value for the particular noise
level. Figure 5.7 illustrates the results of introducing
the second and third noise cases to the measured spectra for
annular flow.

Counts [relative]
Figure 5
.7: Results of introducing the noise and systematic error to the spectrum in
Figure 5.5. a) 10% noise and 25% magnitude decrease; b) 30% noise only.
—1
400
103

104
Table 5.2 summarizes the results of reconstructing the
flow cross section density distribution for the various
combinations of flow regime, void fraction, and noise
content. The reconstructions are rated according to two
criteria: void fraction and image quality. The cross
sectional average void fraction is computed using those
pixels in the image which fall within the pipe boundary.
The image quality is a subjective measure of how well the
reconstructed image matches the true flow distribution.
The error in the measured void fractions is consis¬
tently larger for the low void fraction cases. The detector
response is less sensitive to density variations at low void
fractions than at high void fractions. Therefore, a small
error in the detector response results in a larger error in
the density at low void fractions than at high void
fractions.
Another observation from Table 5.2 is that CST is more
sensitive to systematic error in the spectra than to random
error. The robustness of the reconstruction process to
random error is due to the fact that neighboring channels in
a spectrum are not independent with respect to the indivi¬
dual pixel densities. Due to the finite size of the source
and detector even an abrupt spatial density change results
in smooth spectral features. The best this reconstruction
algorithm can do with a noisy spectrum is to find a density
distribution that results in a spectrum which fits the
average trends and not the individual channels. Systematic

105
TABLE 5.2
Results of reconstruction tests.
Flow
Reqime
Actual
Void
Fraction
Noise Case 1*
Noise Case 2
Noise Case 3
Void
Frac.
Imaqe
Void
Frac.
Imaqe
Void
Frac.
Imaqe
0.30
0.22
9
0.68
0
0.26
9
Bubbly
0.50
0.45
9
0.75
0
0.47
9
0.70
0.71
0
0.85
0
0.72
0
0.30
0.24
9
0.70
0
0.28
9
Strat.
0.50
0.45
9
0.76
0
0.49
9
(sym.)
0.70
0.68
9
0.84
0
0.69
9
0.30
0.34
9
0.72
0
0.34
9
Strat.
0.50
0.48
9
0.77
0
0.52
9
(asym.
0.70
0.68
9
0.84
0
0.70
9
0.30
0.26
9
0.68
0
0.30
9
Annular 0.50
0.46
9
0.73
0
0.47
9
0.70
0.66
9
0.82
0
0.66
9
0.30
0.27
9
0.73
0
0.32
9
Invert
0.50
0.47
9
0.78
0
0.50
9
Annular 0.70
0.67
9
0.84
0
0.69
9
Key to table
• Compares well with the actual distribution.
9 Matches the actual distribution except in a few places.
0 Flow is identifiable but densities are incorrect.
0 Flow is unidentifiable.
* There is some inherent noise in the Monte-Carlo
calculation, but this is less than 6% in all cases.

106
error, unlike random error, affects a spectrum in a more
uniform or smoothly varying manner and does not "average
out" in the reconstruction process. Since systematic error
has a greater impact on the reconstructed density distribu¬
tion it should be minimized.
Resolution
The only available information in CST about the density
distribution comes from the scattered photon spectra. In
order for an object to be resolved it must cause a statisti¬
cally significant change in the spectra by its presence.
This provides a convenient working definition of resolution
and a useful perspective for optimizing system parameters
with respect to resolution. The following discussion
identifies those parameters which affect resolution and
indicates in what manner a parameter can be altered to
improve resolution.
Table 5.3 lists the parameters which affect image reso¬
lution and their ideal values. The ideal values are
classified according to whether the value increases the
spatial information content in the spectra or increases the
number of photons detected thereby reducing the statistical
noise. There exist practical considerations not listed in
the table which force the parameter values in a real system
to deviate from the ideal.
A single detector is only able to give explicit spatial
information in the direction normal to its characteristic
isogonic lines. The energy resolution of the detector; the

107
TABLE 5.3
Some of the parameters which affect image resolution.
Parameter
Ideal value
Detector
size
point
(S)
efficiency
100%
(P)
energy resolution
<<0.001
(S)
number of detector
infinite
(S)(P)
detector-to-
object distance
(S)(P)
Source
size
point
(S)
collimator
perfect absorber
(S)
energy
300 keV to 4 MeV
(P)
source-to-
object distance
as close as possible
(S)(P)
Container (pipe)
thickness
as thin as possible
(P)
. density ,
< density of fluid
(P)
size (diameter)
(P)
S - parameter choice affects the spatial information.
P - parameter choice affects the number of photons detected.

108
size of the detector and source; the source energy; and the
relative locations of the pipe, source, and detector all
combine to define an effective spatial isogonic line
frequency (ESILF). The ESILF is the effective number of
isogonic lines per centimeter measured in the direction
normal to the isogonic lines. This number is not a
constant, but varies over an image. Consider the case of a
point detector and point source where the detector has a 1
KeV energy resolution. The resolution of the detector dic¬
tates an isogonic line spacing no closer than 1 KeV. The 1
KeV energy spacing translates to an actual spatial spacing
determined by the scattering geometry and source energy. If
the point detector is replaced by a detector 1 cm in
diameter then the upper bound on the ESILF is 1 per
centimeter regardless of the detector resolution. This is
due to the scattering angle blur. In general, the size of
the source and detector determine the ESILF in the region
close to the source and the detector energy resolution
determines the ESILF in the regions further away.
Ideally the smallest possible source extension and
detector are used in order to enhance the ESILF. There are
no drawbacks to making a source as small and compact as
possible. The limiting factors are the specific activity of
the source material and the source activity required to meet
the statistical noise and framing rate criteria. A typical
10 Ci cesium-137 source in the shape of a ball can be made
with a diameter of ~6 millimeters.

109
There are definite drawbacks to reducing the size of
the detector, however. Making a detector small adversely
affects its intrinsic efficiency. Also, the smaller a
detector the fewer photons it counts due to field-of-view
effects. The size of the detector and source requires
balancing the competing needs of spatial resolution, time
resolution, and the uncertainty in the density values.
The choice of what type of detector to use is very
restricted. Since all the spatial information of the image
is encoded in the photon energies, detectors with the
highest possible energy resolution are required. Liquid
nitrogen cooled, semiconductor detectors are the only viable
choice with resolutions on the order of 0.002 (i.e. the
full-width-at-half-maximum per photon energy). Another
detector feature to consider is the peak-to-Compton ratio.
Counts in the Compton continuum are unusable due to the
multiply scattered photons in the spectrum. In addition the
tCompton continuum of the high energy photons may interfere
with the low energy end of the singly scattered photon
spectrum. A high peak-to-Compton ratio helps to avoid these
problems.
Image Contrast
The discussion on resolution has concentrated on what
might be called the spatial sampling rate. An object
smaller than the spatial sampling rate may enter the image
plane and still cause a statistically significant change in

110
the spectra. The object, however, is blurred in the recon¬
structed image. It is also possible for an object to be
larger than the spatial sampling rate and not cause a
statistically significant change in the spectra. This can
occur for two reasons. First, the object could have the
same electron density as the rest of the surrounding
material and be indistinguishable from it. Second, if the
object is sufficiently shielded by the bulk material, then
the detectors never see the scattered photons coming from
the region of the object. The relationship between the size
and density of an object which can be sensed is analogous to
the relationship between sharpness and contrast in
photography.
The two factors affecting whether an object can be seen
are its size and the difference between its electron density
and that of the bulk material. Specifying the minimum size
and density disturbance that can be sensed necessitates
estimating what constitutes a statistically significant
change in a spectrum. The two major sources of error which
define statistical significance are measurement-data
processing error and reconstruction error. The measurement-
data processing error has three main components: counting
statistics, background removal, and detector efficiency
correction. Equation (2.8) gives an expression from which
the error in these three components can be related to the
error in the singly scattered photon spectrum. The relative

Ill
error, xcr in the single scatter counts for a particular
channel, c, in a spectrum is given by
J a
- -Í -t
'D
1
o
+
Si,2 a2
1/2
] /
(5.1)
where singly scattered photon spectrum,
a£ is the standard deviation of the detector
efficiency correction,
aD is the standard deviation of the detected count in
channel c, and
aB is the standard deviation of the background count
in channel c.
The variance of the total detected counts, Dc, and the
background counts, Bc, can be approximated by
= Dc and aB = Bc . (5.2)
The efficiency correction, e, does not represent a signif¬
icant source of error. This is because the value of the
efficiency correction is not as important as its functional
form. Any error in the value of the efficiency correction
can be compensated for by a calibration procedure, but only
as long as the correction in one channel relative to the
others is done correctly. The functional variation of
detector efficiency with energy is well understood in the
range of energies of interest in Compton scattering. The
relative error in the efficiency correction is probably on

112
the order of 1%
i.e. a£/e =0.01. (5.3)
The reconstruction error involves the error in the
cross section library, the usual error in finite precision
arithmetic, and the error associated in representing the
two-phase flow with discrete pixels. Hubbell [34] estimates
the error in the National Bureau of Standards published
cross section library for the energy range of interest (200
keV to 4 MeV) to be on the order of 1%. The reconstruction
algorithms developed here are not very sensitive to round
off error, etc. and with the advent of microcomputers which
use 64 bit real arithmetic means that this error can be
neglected. The primary source of error in the reconstruc¬
tion process is in the discrete representation of the flow.
Experience has shown that most of this error is in calcula¬
ting the attenuation factors. Using small pixels (~0.1
mean-free-path) and averaging the attenuation factors from
several random points per pixel can reduce this error to the
order of 5%. This error mainly affects those pixels on the
edge of features in the flow. The overall result of the
reconstruction error is that even when the correct density
distribution (in discrete form) is used the difference
between the calculated spectra and the measured spectra is
about 5%
aF/sl =0.05 .
i. e.
(5.4)

113
The change in the detected number of photons, ADC,
which an object must cause in order to be seen is given
approximately by
ADc/ec > {a g + Qp}1/2 (5.5)
2
where ap is the variance of the reconstruction error. The
value of ADC for a particular object can be estimated using
a modified form of equation (2.1)
ADC = S°ApywAzAiA0ecPSFc (5.6)
where Ap is the difference between the object density and
the bulk material density and Az is the mean distance
through the object as seen by the source photons. The
attenuation function, A^, and the isogonic segment charac¬
teristic function, PSFC, are both functions of the object
size, shape, and orientation. Even assuming a simple object
shape, PSFC is not easily expressed as an analytic function
of Az and, in general, PSFC must be calculated numerically.
In practice, equation (5.6) is computed for several sizes of
an object and the value of ADC which satisfies equation
(5.5) is used to interpolate the results to find ApAz.
The simplicity of equations (5.1) through (5.6) belie
the problem of specifying a minimum size object and density
variation that can be sensed. The attenuation factors in
equation (5.4) introduce spatial, void fraction, and flow

114
regime dependence to ApAz. The detector efficiency adds an
energy dependence and the other variables introduce yet
other functional dependencies.
Note that the choice of source energy is coupled to the
problem of image contrast by the mass attenuation and scat¬
tering coefficients. The mass interaction coefficients
decrease as the source energy increases. At higher source
energies the photons penetrate better and shielding by the
bulk material is less of a problem. At lower source
energies the photons interact better and there is better
sensitivity to a given change in the electron density.
There exists, for each particular imaging system, a source
energy which optimizes the detector response by balancing
the effects of attenuation and interaction. An optimum
source energy can be found which maximizes the number of
detected photons coming from a point in the flow. That is,
a search is made for the source energy which satisfies
3D(x,y)/3E = 0 . (5.7)
In the case of the experimental setup in Figure 5.3 the
optimum source energy is close to that of cesium-137 for 0%
void fraction and a point in the flow furthest from the
source and detector. The optimum energy decreases with
increasing void fraction and/or for points closer to the
source and detector.

115
The choice of source energy must also take into account
the practical considerations of finding an isotope with a
mono-energetic gamma of the right energy and with a
reasonable half-life. A dual energy source such as cobalt-
60 could be used with some care to separate the two singly
scattered photon spectra. Another practical problem is
source handling and shielding. High energy, high activity
sources are impractical due to the amount of shielding
material required.
Image contrast is also affected by the presence of a
dense scattering and attenuating object such as a steel
pipe. The pipe wall introduces a threshold effect. If the
pipe is thick enough and its electron density much greater
than the bulk material being imaged, then the scattered
spectrum of the pipe wall can dominate the detector
response. The statistical variation in the pipe wall
spectral component sets a threshold of statistical signifi¬
cance which variations in the bulk material must overcome in
order to be sensed. The counts in the singly scattered
photon spectrum due to the bulk material, S^-^, must be
greater than the expected Poisson statistical variation in
the spectrum,
i-e- Sbulk > aS (5*8)
where ag is given by equation (5.1). The threshold effect
is important at high void fractions where there are few

116
photons interacting in the bulk fluid compared to the pipe.
Increasing the data acquisition time in order to increase
the number of detected photons from the bulk fluid is one
way of reducing the impact of the threshold.
Figure 5.8 shows the spectrum from a 50% void fraction
annular flow (same as Figure 5.5) with a 0.6 centimeter
thick iron pipe taken into account. The peaks at the upper
and lower ends of the single scatter spectrum are dominated
by the pipe wall and it is more difficult to see changes in
the bulk material density in the isogonic segments
associated with these two peaks. The pipe also attenuates
the singly scattered photons from the bulk fluid especially
toward the low energy end of the spectrum.
A CST instrument is a highly integrated system in which
each design parameter is chosen in concert with the others
and with respect to the characteristics of the object being
imaged. The design process can be understood as optimizing
three basic system specifications: accuracy, image resolu¬
tion, and contrast resolution. Accuracy determines how
close a reconstructed image matches the gross features of
the object being imaged. Image resolution and contrast
resolution determine the minimum size features of an
object which can accurately be reconstructed in the
image.

117
Energy [keV]
Figure 5.8: Detected spectrum from a 50% void fraction
annular flow inside a 0.6 cm thick walled
steel pipe.

CHAPTER 6
SUMMARY AND CONCLUSIONS
This study investigates and demonstrates several
methods of inferring the flow parameters of two-phase flows
using wide beam illumination and few (i.e. 2) detectors.
The first step is the recognition of the singly scattered
photon spectra as projections of the phase distribution.
The spatial distribution of the fluid-vapor phase in a slice
of the pipe is encoded with respect to energy in the singly
scattered photon flux. Two basic methods are detailed
for decoding the spatial information: the method of
spectral moments and the method of computed tomography.
Previous authors have recognized that tomography of '
two-phase flow is possible using Compton scattering of gamma
rays. The majority of their efforts employ elaborate colli-
mation schemes to reduce the likelihood of multiply
scattered photons reaching the detectors. This report
demonstrates using 3-D Monte-Carlo photon transport calcula¬
tions that the multiply scattered photon spectrum for a wide
beam source and uncollimated detector setup lacks any
structure which would interfere with extracting the singly
scattered photon spectrum. A simple background subtraction
procedure appears sufficient for removing the multiply
scattered photon spectrum component. Not only is the
118

119
detector collimator unnecessary it is undesirable since it
can reduce the sensing of singly scattered photons and
thereby degrade system performance.
Summary
Examination of the low-order moments of the singly
scattered photon spectra for a one source, two detector
configuration provides sufficient information for flow
regime identification. The two detectors are placed in
symmetric positions with respect to the source and pipe.
Flow asymmetries are revealed by differences in their
spectral moments. Further identification of the actual
symmetric or asymmetric flow pattern is made on the basis of
the values of the first and second moments of the spectra.
Quantitative measures of the flow may be possible using this
method and a spectral calibration technique.
The real'focus of this report is the adaptation and
successful demonstration of computed tomography techniques
with Compton scattering. Three series expansion techniques
are chosen for adaptation because of their iterative
nature— ART, SIRT, and ILS. The adaptation process
involves deriving the equations for each algorithm from the
forward Compton scattering relations and altering the
structure of the iterations to include a recalculation of
the transport matrix. Recalculation of the transport
matrix overcomes the nonlinearity associated with
attenuation without having to include it explicitly in any
derivatives. Application of the modified ART, SIRT, and ILS

120
algorithms to a hot-spot and a cold-spot reconstruction
problem indicates that SIRT and ILS are more accurate than
ART. A more extensive testing of the ILS algorithm for a
variety of model flow regimes demonstrates the potential of
CST as a quantitative measurement technique.
Monte-Carlo backprojection is another reconstruction
technique looked at in this report. This is a new recon¬
struction method which operates by statistically cross
correlating the spectra of a few detectors to produce a
tomographic image. This technique appears well suited to
situations where the spectra contain a lot of structure.
One example would be the binary type flow patterns expected
for high frame rate imaging of two-phase flows. Another
application would be to identify and localize flow anomalies
by looking at the difference between a time averaged
spectrum and an instantaneous spectrum. The difference
spectrum would show structure due to any flow anomalies.
Based on the observations and results of Chapter 5 a
CST system is best designed with specific measurement goals
in mind. The design process involves striking a balance
between image resolution, image accuracy, framing rate, and
sensitivity. The balancing or optimization is such that
there'is a tight coupling between a CST system and a
particular application. Even with careful optimization
it may not be possible to achieve set design goals due to
practical limitations.

121
A Practical Application
As a practical illustration and summary of the design
concerns discussed in Chapter 5 consider the problem of
imaging the phase distribution of a steam-water flow inside
a 4 inch schedule 40 steel pipe. The source-pipe-detector
configuration is assumed to be the same as shown in Figure
5.3. Figures 6.1 and 6.2 show the results of a Monte-Carlo
calculation of the scattered photon spectrum for this setup
when the void fraction is zero. The two quantities of
interest to an experimenter are the source activity required
to achieve a set framing rate and the minimum size bubble
that can be sensed.
Estimating the required source strength necessitates
first calculating the minimum number of photons which must
be counted in a single channel in order to overcome
statistical noise in the measurements. The source strength
is given by dividing the minimum number of photons for a
single channel, Dc, by the number of photons detected in
channel c per emitted source photon and multiplying by the
framing rate. A relation for Dc can be obtained by
substituting equations (5.1) and (5.2) into equation (5.8)
rD(SNR)2(rD + rfi)
Dc = 5 5—
e(l - r2(ae/e)2)
where rD is the ratio Dc/S¿,
rB is the ratio Bc/S¿,
rc is the ratio S¿/S¿,
(6.1)

Photons detected per source photon [photons/Ci-sec
Figure 6.1: Results of a Monte-Carlo calculation of the scattered photon spectrum for a
4 inch schedule 40 pipe full of water. The calculation assumed a 100%
efficient detector and a 1 keV MCA channel width setting.
122

123
Energy [keV]
Figure 6.2: Detail of the enclosed region on Figure 6.1
showing the total number of photons detected, D ,
the multiply scattered photons, Bc, the singly c
scattered photons, S¿, and the singly scattered
photons from the bulk fluid, S¿.

124
SNR is the signal-to-noise ratio, and
Dc, S¿, Bc, and are defined pictorially in Figure 6.2.
Equation (5.8) gives an inequality which must be satisfied
in order to extract information about the flow from the
photon spectra. Replacing the inequality by an equality is
equivalent to specifying a working SNR of 1. Figure 6.3
illustrates some flow reconstructions obtained when the SNR
is near 1.
The ratios rD, rB, and rc for channel c estimated from
Figure 6.2 are 2.79, 1.27, and 1.51, respectively. Substi¬
tuting these ratio values, equation (5.3), an assumed
detector efficiency of 15%, and an assumed SNR of 2 into
equation (6.1) indicates that 300 photons must be counted in
channel c (i.e. Dc = 300) per measurement. The Monte-Carlo
calculations indicate that 36 photons are detected in
channel c per Curie of source activity per second for a 100%
efficient detector. The required source activity for a 15%
efficient detector is given by
(300 photons/frame)(10 frames/sec)
activity = = 556 Ci
(0.15) (36 photons/Ci-sec) (6.2)
The minimum size steam bubble that can be sensed with
this CST system is found using equations (5.5) and (5.6).
The RHS of equation (5.5) is evaluated using equations (5.1)
and (5.4) along with the information from the source
activity calculation which yields

125
Figure 6.3: Reconstruction of a 50% void fraction annular
flow regime for three different SNR values of
photon spectra.

126
ADC > 23 counts . (6.3)
The same comment about the SNR with respect to equation
(6.1) applies to this inequality and it is assumed that ADC
needs to be twice as large as the RHS of equation (6.3)
i.e. ADC = 46 or ADC/DC = 46/300 = 0.15 . (6.4)
Figure 6.4 plots the results of a numerical evaluation of
equation (5.6) for spherical bubbles located in the center
of the 4 inch pipe. The peak in the curve arises because
of the decreased attenuation of the scattered photons as
they pass through the larger bubbles. The smaller bubbles
do not significantly affect attenuation and simply result
in less photons being scattered. A conservative estimate
of the minimum size steam bubble is read from the graph as
3.5 cm. In summary a CST system to measure steam-water
flow in a 4 inch schedule 40 pipe at 10 frames per second
requires a 556 Curie cesium-137 source and is able to
sense bubbles on the order of 3 cm in diameter.
Conclusions and Future Research
Looking specifically at the two-phase flow measurement
problem there are several conclusions which can be drawn
from the analysis of this report. It is possible to build a
CST system to image two-phase flow inside steel pipes using
a single source and two detectors. Such a CST system can be
expected to function well in the presence of moderate

ADC/D
Figure 6.4: Results from the numerical evaluation of equation (5.6).
127

128
statistical noise, but may give erroneous results if system¬
atic errors exist in the data. The performance of the
system is best at low void fractions and degrades with
increasing void fraction. The pipe can place a practical
upper limit on how high the void fraction can go and
meaningful measurements still be made. If the pipe wall is
too thick it may not be possible to design a functioning
system at all. Some of these limitations can be over come
by augmenting the CST design of this report with a few
detectors to measure the attenuation of the transmitted
source beam. The attenuation measurements complement the
scatter measurements since they perform best at high void
fractions. The potential advantages of a hybrid scatter-
transmission system warrants further investigation and it is
concluded that the next logical step would be the construc¬
tion of a hybrid CST prototype system. A reasonable design
goal is to image air/water of steam/water flows inside thin
walled pipes with diameters in the range of 8 to 16 cm.
Real-time data analysis and display could be achieved by
integrating a super-micro computer into the system. Besides
functioning as a useful laboratory instrument the prototype
would provide a means for evaluating analytic techniques for
boosting the effective system framing rate without requiring
excessively large source strengths. Despite the evidence in
this report that CST is a viable imaging technique, a
prototype system is necessary to gain the acceptance of
other researchers, particularly in the medical community.

129
The basic CST concepts discussed here require mono-
energetic sources of gamma rays. This limits the choice of
sources to only a few commercially available isotopes.
Extension of the techniques to poly-energetic sources would
allow closer optimization with respect to source energy and
perhaps the use of X-ray machines. The higher photon
outputs of X-ray machines appear necessary to achieve true
stop action imaging.

APPENDIX I
MONTE-CARLO TRANSPORT PROGRAM
The 3-D Monte-Carlo photon transport computer code is
written in Microsoft FORTRAN for the IBM PC. This program
is the one used to generate the "measured" data of this
report. The program is written for a very specific
geometry: a cylindrical pipe running parallel to the z-axis
and illuminated by an external source located at the origin.
The histories of a fixed number of photons are followed
from birth through a fixed number of scattering
interactions. It is assumed that only Compton scattering
collisions and photoelectric absorption interactions take
place. Three different source models are available: point,
distributed (i.e. cylindrical), and pencil beam. All
photons are born at the origin with the point source and
pencil beam models. The distributed source model uniformly
chooses the coordinates of a photon's birth from within a
cylindrical source. The initial trajectory of all photons is
chosen so that their paths cross the pipe. Any trajectories
allowed by the source collimator which do not intersect the
pipe are ignored. The point of intersection between a
photon's initial trajectory and the pipe wall is used to
start the photon tracking process.
130

131
The photon tracking process involves the repeated exe¬
cution of three steps: 1) choose an interaction point
inside the pipe along the current photon trajectory, 2)
score the fraction of photons interacting at this point
that would reach the detector, and 3) choose a new photon
trajectory by sampling the Klein-Nishina differential
scattering cross section. These three steps are executed
once for each order of scattering. The importance of each
photon is corrected in step 1 for the probability of the
photon not escaping the pipe and for the probability that
the photon is scattered instead of absorbed. The correction
for the probability of not escaping is necessary since all
photons are forced to have their next interaction inside the
pipe. The scoring in step 2 takes into account the photon
importance.
It is often stated that Monte-Carlo does nothing more
than trade the difficulties of solving a transport equation
for an equally difficult geometry problem. Figure 1.1
describes a common geometry problem which must be solved by
the program and should be helpful in understanding some of
the logic and equations in the subroutines. An effort has
been made to use variable names consistent with this figure.
Two floppy disks containing this computer code are
available by contacting the author through the Department of
Nuclear Engineering Sciences, 202 Nuclear Sciences Center,
Gainesville, Florida 32611.

= R - (R sin9)
P
COS0 = V *v
P o
Figure I .1: Geometry and variable definitions used in the transport program.
132

APPENDIX II
INSIGHTS AND FINER POINTS
Parallel Versus Diverging Source Beams
The reconstruction algorithms in Chapter 4 explicitly
ignore the density terms in the arguments of the exponential
attenuation factors. The expression for the number of
singly scattered photons counted in a single MCA channel can
be simplified to the following form
Sc = pc*Constant*exp[ -p'*(uiLi + yQL0)] (II.1)
where pc is the average density of the isogonic segment
associated with channel c,
p' is the average density in the flow cross section
as seen by the incoming source and outgoing
scattered photons,
Pj_,U0 are the mass attenuation coefficients for the
incoming source and outgoing scattered photons,
respectively, and
Lj_,LQ are the average distances traveled in the two-
phase flow by the incoming and outgoing photons,
respectively.
In a fan beam source geometry pc and p' are not necessarily
related since the paths of the incoming source photons and
the outgoing scattered photons only traverse the isogonic
segment c for relatively short distances (see Figure 2.5).
The reconstruction algorithms take advantage of the
133

134
observation that when an individual pixel density is
adjusted the associated attenuation factor for that pixel
does not change. Obviously the attenuation factors for the
other pixels do change if the incoming or scattered photons
of those pixels cross the pixel whose density is changed.
In the case of a parallel source beam the path of an
outgoing scattered photon is colinear with the isogonic
segment (see Figure 2.6). Equation (II.1),' rewritten for a
parallel source beam, is
Sc = Pc*Constant*exp[ -(p'y^L.^ + PCU0L0)] . (II.2)
A reconstruction algorithm can more easily account for the
linear-exponential dependence on pc in equation (II.2) than
in equation (II.1).
Figures II.1 and II.2 graphically depict equations (II.1)
and (II.2) for an isogonic segment crossing the middle of the
pipe in Figure 5.3. Equation (II.2) can exhibit a maximum
with respect to pc if the pipe is large enough. The
presence of a maximum can lead to an ambiguity in
determining pc. The upshot of this is that the diverging
source beam configuration leads to a more tightly coupled
set of reconstruction equations than the parallel source
beam.

Figure II.1: Detector response for channel c versus the isogonic segment average
density, p , and the average density seen by the source and scattered
photons, p". Wide source beam illumination is used.
135

Figure II.2:
Detector response for channel c versus the isogonic segment average
density, p , and the average density seen by the source photons, p'
Parallel sSurce beam illumination is used.
136

137
Multiple Scatter Versus Detector Position
A convenient measure of the importance of the multiply
scattered photons is the ratio of the singly scattered
photons to multiply scattered photons averaged over the
singly scattered photon spectrum. This ratio varies with
detector position for a lot of reasons. The most important
reason is that Compton scattering is non-isotropic. Figure
II.3 shows the same experimental setup as Figure 5.3, but
with three different detector positions. Figure II.4 shows
the detected spectra for those positions. The single-to-
multiple scatter ratio is highest for position 3 and lowest
for position 2.
Position 3 is favored for reducing the importance of
multiple scattering, however, it may not be the best choice.
Using position 3 would require a higher energy resolution
detector in order to have the same ESILF as position 2 or 3,
since the single scatter spectrum is narrower in position 3.
Reference [35] discusses similar ways of optimizing the
detector placement.
Other factors to consider in placing the detector are
detector efficiency versus angle of incidence, geometric
attenuation, and spatial resolution. It is desirable to
place the detector as close as possible to the object in
order to increase the number of detected photons and the
width (on the energy axis) of the singly scattered photon
spectrum. The wider the spectrum the more spatial
information is available for a given energy resolution

(12,15)
D
(23,13)
y
A
x, y)
(4,9)
Figure II.3: Three detector positions used to illustrate the variation of the multiple
scatter component of the spectrum with detector position.
138

Counts [relative]
Figure II.4:
Spectra for the detectors in Figure II.3.
139

Figure II.4—continued.
Detector 2
Single-to-multiple scatter ratio
7.67
Energy [keV]
600
140

Counts [relative]
Figure II.4--continued.
141

142
detector assuming the detector and source sizes are
negligible.
The detector efficiency can be a strong function of a
photon's angle of incidence to the detector. This is a
problem when photons of the same energy (i.e. coming from
the same isogonic segment) have widely ranging angles of
incidence. The usual efficiency correction procedures
assume all photons enter at the same angle and can not
correct for such a variation. Moving the detector away from
the object reduces the range of angle of incidence. Planar
detectors are more sensitive to angle of incidence and less
efficient than coaxial detectors. Therefore, it is recom¬
mended that coaxial detectors be used in CST.
Source Placement
There are some practical matters to consider in the
placement of the source. Photon economy is directly related
to the source-to-object distance. The closer the source is
to the object the larger the ratio of photons illuminating
the object to the photons emitted. The trade-off in moving
the source closer to the object is some loss in flexibility
in the placement of the detectors. The detectors must be
placed (or shielded) such that they do not have a line-of-
sight view of the source.
Another point to consider is backscatter from any sur¬
rounding structures. The wide source beam associated with
small source-to-object distances is more likely to cause

143
problems by illuminating nearby objects. Scatter from these
nearby objects could interfere with the measurements.
Convergence Criterion
The choice of the criterion used to end the reconstruc¬
tion iterations can be tricky. Reference [36] looks at a
similar situation where there are only a few projections
available and considers two types of criteria:
a) criteria based on properties of the reconstruction,
and
b) criteria based on the residual errors in the
projections.
Their conclusion based on tests is that type b) methods are
better when there are only a few views.
The type b) convergence criteria for CST corresponds to
examining the difference between the measured spectra and
the calculated spectra for the reconstruction. Appendix V
contains' plots of t'he spectrum difference function and the
discrepancy function (see the discussion of Figure 4.11)
versus iteration number. These plots show that there is a
very weak correlation between the rate of change in the
spectrum difference function and how close the reconstruc¬
tion is to the actual flow pattern as measured by the dis¬
crepancy function. This author prefers to use a fixed
number of iterations instead of a convergence criterion.

APPENDIX III
RECONSTRUCTION COMPUTER PROGRAMS
The reconstruction algorithms derived in Chapter 4
are implemented here in Turbo Pascal for an IBM PC.
There is one main program for each of the four
algorithms: Monte-Carlo Backprojection (MCB), ART, SIRT,
and ILS. There are also four separate procedures which
are the actual implementation of these algorithms. All
remaining procedures used are common to the four main
programs and are included at compile time.
Running the Programs
Prior to running the programs for a particular source-
pipe-detector configuration a Pixel Spread Function (PSF)
file must be generated for that particular setup. The PSF
file contains the values of the PSF matrix introduced in
Chapter 4. The program Make_PSF is used to generate this
matrix using a Monte-Carlo technique and to store it in the
proper format.
The programs are designed to be run interactively. The
first thing the programs do is read in the input data file
and the PSF file. Reading in the PSF file takes approxi¬
mately 3 minutes. The programs then run through 10
144

145
iterations of their respective algorithms (~30 minutes).
The programs then prompt for the name of an output file.
The input data file contains the information about the
system geometry, descriptions of the source and detector,
and the measured spectrum for each detector. The data
fields are free format, however, data types must be adhered
to. The file SAMPLE.INP contains enough descriptive
comments that it can be used as a model for other problems.
Three floppy disks are available containing these
programs plus sample data by contacting the author through
the Department of Nuclear Engineering Sciences, 202 Nuclear
Sciences Center, Gainesville, Florida 32611.
Flow of Programs
The ART, SIRT, and ILS programs all have the same basic
internal structure. The first statement reserves space on
the heap for the PSF function. The next block of statements
initialize all the global variables, loop counters, etc.
Following the initialization code block is a repeat-until
block which executes the particular reconstruction algorithm
a fixed number of iterations. The final statements save the
reconstructed density distribution and other results in an
output file.

APPENDIX IV
MAXIMUM ENTROPY METHOD
The maximum entropy method begins by redefining the
fluid density distribution as a probability density
distribution. The density distribution, p, usually sought
is replaced by the probability function, E, where
it k = i
and (IV.l)
fk i 0.
The relationship between f and p is
fk = Pk/(P*K) (IV.2)
where p is the average density and K is the total number of
pixels. The average density may not be known exactly and
changing the normalization of f to an inequality allows
using the fluid density in place of the average density.
The entropy of a probability distribution is
H(f) = -Ifkln(fkK) . (IV.3)
The derivation given here is based on the work of
Minerbo [28].
146

147
The maximization is carried out by introducing a Lagrange
multiplier XCT for each constraint condition (i.e. the
y
measured data) and forming the Lagrangian
'F(f) = H (f ) - £xg* (Sg - s| ( f ) ) (IV.4)
Taking the functional derivative with respect to f while
assuming that the transport matrix [equation (4.1)] is
independent of f and setting it equal to zero yields
H-l - ln(fkK) + (l/PfiuidK)£TgkV = ° • (IV*5)
k g
Now it is argued that in order to have K linearly indepen¬
dent density values requires that the term in brackets be
zero for each value of k. Solving for fk yields
fk = exp{ 1 - UgTgk}/K . (IV. 6)
g
The Xg's are found by substituting the density values given
by equations (IV.6) and (IV.2) into the constraint conditions
[i.e. equation (4.1)]. This gives G equations of the form
Sg ' = Pfluid* £Tgk*exP{1 “ IXgTgk}
k g
where the - values are the measured data.
y
(IV.7)

148
It is interesting to note here the relationship between
the maximum entropy solution and the least squares solution
[see equation (4.28)]. When the number of pixels, K, is
less than the number of measured data points, G, the lambdas
are undefined. When the number of pixels is the same as the
number of data points the maximum entropy solution is iden¬
tical to the least squares solution with the weights set
equal to one.
When the number of pixels is larger than the number of
data points the maximum entropy solution exists. The infor¬
mation needed to define the extra pixels is provided by the
entropy condition. The problem is that a nonlinear set of
equations must be solved for the lambdas. The experience of
this author with using the maximum entropy method under
these conditions is that the solution fails to converge.
Even for K only slightly greater than G leads to an ill-
conditioned set of equations. The problem seems to be that
the f^ values are too sensitive to small changes in the
lambdas coupled with the fact that the transport matrix
elements are functions of the density.
An even more rigorous derivation of the maximum
entropy solution can be made by explicitly including the
dependence of the transport matrix on the density distribu¬
tion. The equation analogous to equation (IV.5) arrived at
by this route leads to a nonlinear relation between the f's
and lambdas which just complicates the solution even more.

APPENDIX V
DATA
This appendix is divided into two parts. The first
part contains a representative selection of the data used in
compiling Table 5.2. Figures V.l through V.15 are some of
the spectra generated by the Monte-Carlo transport code.
Figures V.16 through V.29 are plots of the reconstructed
flow distributions for various cases of flow regime, void
fraction, and noise content. Figures V.30 through V.34 are
plots of the spectrum difference and discrepancy functions
showing the different types of convergence behavior.
The second part contains plots of the spectra from
actual laboratory experiments. This data is not usable
because of a lack of available sources to calibrate the
detector and electronics. Figure V.35 shows a sketch of the
experimental setup used. Figure V.36 illustrates the steps
of correcting the spectrum for the detector efficiency and
removing the multiply scattered photon spectrum. The
multiple-scatter "background" is removed by fitting a 2nd
order equation to the two small spectrum segments on either
side of the single scatter peak.
Also, shown in Figure V.37 is the spectrum generated by
the Monte-Carlo photon transport code for the same annular
setup. The slight misalignment of the two spectra is
149

150
probably due to poor energy calibration in the MCA and some
difficulty in determining the detector efficiency. A three
point energy calibration would be better than the two point
procedure actually used. The detector efficiency was
determined using sources of unknown pedigree. Although
energy calibration of an MCA and detector efficiency deter¬
mination are routine procedures in most counting laborato¬
ries this figure illustrates the need for care in performing
these tasks.

Counts [relative]
<1)
>
•H
-p
nJ
r—I
Q)
U
W
-P
a
d
o
u
Figure V.l:
Detected spectra from a 30% void fraction annular flow regime. a) Spectrum
of upper detector; b) Spectrum of lower detector.
151

Counts [relative]
Figure V.2: Detected spectra from a 50% void fraction annular flow regime. a) Spectrum
of upper detector; b) Spectrum of lower detector.
152

Counts [relative]
Figure V.3: Detected spectra from 70% void fraction annular flow regime. a) Spectrum
from the upper detector; b) Spectrum from the lower detector.
153

Counts [relative]
Figure V.4:
Detected spectra from a 30% void fraction inverted annular flow regime,
a) Spectrum from the upper detector; b) Spectrum from the lower detector.
154

Counts [relative]
0.1
a)
>
•H
-P
m
i—i
0
w
4J
£
£
O
O
200 Energy [keV]
(b)
"1
400
Figure V.5:
Detected spectra from a 50% void fraction inverted annular flow regime,
a) Spectrum of upper detector; b) Spectrum of lower detector.
155

Counts [relative]
1
400
Figure V.6:
Detected spectra from a 70% void fraction inverted annular flow regime,
a) Spectrum of upper detector; b) Spectrum of lower detector.
156

Counts [relative]
0.1
(U
>
•H
-P
rtf
r—4
o
W
-P
C
3
O
U
(b)
“1
400
Figure V. 7 :
Detected spectra from a 30% void fraction bubbly flow regime. a) Spectrum
of upper detector; b) Spectrum of lower detector.
157

Figure V.8: Detected spectra from a 50% void fraction bubbly flow regime- a) Spectrum
of upper detector; b) Spectrum of lower detector.
400
158

Counts [relative]
0
Figure
.9: Detected spectra from a 70% void fraction bubbly flow regime. a) Spectrum
of upper detector; b) Spectrum of lower detector.
159

Counts [relative]
0.1
CD
>
+j
Ctf
r—I
CD
M
W
4J
G
G
O
U
0
(b)
Figure V.10: Detected spectra from a 30% void fraction stratified symmetric flow regime,
a) Spectrum from upper detector; b) Spectrum from lower detector.
400
160

Counts [relative]
>
-p
(0
i—i
0)
u
m
-p
c
3
o
o
Figure V.11:
Detected spectra from a 50% void fraction stratified symmetric flow regime,
a) Spectrum from upper detector; b) Spectrum from lower detector.
161

Counts [relative]
Figure V.12: Detected spectra from a 70% void fraction stratified symmetric flow regime,
a) Spectrum from upper detector; b) Spectrum from lower detector.
162

Counts [relative]
Figure V.13: Detected spectra from a 30% void
a) Spectrum from upper detector;
fraction stratified asymmetric flow regime,
b) Spectrum from lower detector.
163

Counts [relative]
Figure V.14: Detected spectra from a 50% void fraction stratified asymmetric flow regime,
a) Spectrum of upper detector; b) Spectrum of lower detector.
164

Counts [relative]
Figure V.15:
Detected spectra from a 70% void fraction stratified asymmetric flow regime,
a) Spectrum of upper detector; b) Spectrum of lower detector.
165

166
Annular flow regine as it should appear when
correctly reconstructed, a) 30% void
fraction; b) 50% void fraction; c) 70% void
fraction.
Figure V.16:

167
Figure V.17: Reconstruction of annular flow for noise case 1.
a) 30% void fraction; b) 50% void fraction;
c) 70% void fraction.

168
Figure V.18: Reconstruction of annular flow for noise case 2.
a) 30% void fraction; b) 50% void fraction;
c) 70% void fraction.

169
Figure V.19: Reconstruction of annular flow for noise case 3.
a) 30% void fraction; b) 50% void fraction;
c) 70% void fraction.

170
Inverted annular flow regime as it
appear when correctly reconstructed
void fraction; b) 50% void fraction
void fraction.
should
. a) 30%
; c ) 7 0 %
Figure V.20:

171
(c)
Figure V.21: Reconstruction of the inverted annular flow for
noise case 1. a) 30% void fraction; b) 50%
void fraction; c) 70% void fraction.

172
Figure V.22: Bubbly flow regime as it should appear when
correctly reconstructed, a) 30% void fraction;
b) 50% void fraction; c) 70% void fraction.

173
Figure V.23: Reconstruction of the bubbly flow for noise
case 1. a) 30% void fraction; b) 50% void
fraction; c) 70 % void fraction.

174
(c)
.BBT"'
Jig*
ty'.:;'.
•?Aíf.
V-1
UM.W
•;H{
M?S_. . ....
. ¡?T;r, í/ütVfVÁ ;•:
\y-C...«5 ,:..\/ •.,
.ÍXÍ* t,’.
•'■?'• ví •.
v«í
‘ Vy-CS-A
*S T.W>\>.
Figure V.24: Reconstruction of the bubbly flow for noise
case 2. a) 30% void fraction; b) 50% void
fraction; c) 70% void fraction.

175
Figure V.25: Reconstruction of the bubbly flow for noise
case 3. a) 30% void fraction; b) 50% void
fraction; c) 70% void fraction.

176
Stratified symmetric flow regime as it should
appear when correctly reconstructed, a) 30%
void fraction; b) 50% void fraction; c) 70%
void fraction.
Figure V.26:

177
(b)
(c)
Figure V.27:
Reconstruction of the stratified symmetri
for noise case 1. a) 30% void fraction;
void fraction; c) 70% void fraction.
c flow
b) 50%

178
Stratified asymmetric flow regime as i
appear when correctly reconstructed,
void fraction; b) 50% void fraction; c
void fraction.
t should
a) 30%
) 70%
Figure V.28:

179
(b)
(c)
Figure V.29:
Reconstruction of
for noise case 1.
void fraction; c)
the stratified asymmetric flow
a) 30% void fraction; b) 50%
70% void fraction.

Figure V.30: Example of smooth convergence.
180

1
cu
o
tí
cu
tí
O
4-4
4-4
•r-í
Q
e
d
tí
-p
o
(U
Ch
cn
0
Iteration
“I
10
Figure V.31: Example of overiteration with resulting slight divergence.
181

Discrepancy
Figure V.32:
Example of overiteration with resulting divergence.
182

Figure V.33: Example of overiteration with oscilating divergence.
183

Figure V.34: Example of failure to converge.
184

Coax
Detector
Figure V.35: Experimental source-pipe-detector configuration.
185

6-
5-
(a)
Figure V.36: Experimental Spectrum from the annular plastic phantom. a) raw data;
b) detector efficiency corrected data with dashed line indicating fit
to multiple-scatter background.
186

Figure V.36—continued.
1 i 1 1 "I
400 600
Energy [KeV]
(b)
187

Figure V.37: Comparison of Monte-Carlo generated spectrum to experimental spectrum from
the annular flow phantom.
188

7
Figure V.38: Experimental spectrum from a stratified flow phantom.
189

Counts
7
Energy [KeV]
Figure V.39: Experimental spectrum from a bubbly flow phantom.
190

REFERENCES
1. Hewitt, G.F., Measurement of Two-Phase Flow Parameters,
Academic Press, London (1978).
2. Heiselmann, H.W., Davis, L.L., "Measurement Techniques
for Transient Two Phase Coolant Flow During Blowdown,"
Proceedings of Two-Phase Flow and Heat Transfer
Symposium-Workshop, Fort Lauderdale, Florida, (October
18-20, 1976), (distributed by Pergamon Press).
3. Delhaye, J.M., Cognet, G., "Measuring Techniques in Gas-
Liquid Two-Phase Flows," Proceedings of International
Union of Theoretical and Applied Mechanics Symposium,
Nancy, France, July 5-8, 1983, Springer-Verlag, Berlin
(1984).
4. Bergles, A.E., Collier, J.G., Delhaye, J.M., Hewitt,
G.F., Mayinger, F., Two-Phase Flow and Heat Transfer in
the Power and Process Industries, Hemisphere Publishing
Company, Washington, D.C. (1981).
5. Bryant, D.T., Payne, J.A., Firmin, D.N., Longmore, D.B.,
"Measurement of Flow with NMR Imaging Using a Gradient
Pulse and Phase Difference Technique," Journal of
Computer Assisted Tomography, Vol. 8, pg. 588 (1984).
6. Jacobs, A. M., "Nonintrusive transaxial tomography
technique for velocity profile measurement," Proceedings
of the Society of Photo-Optical Instrumentation
Engineers, Vol. 182 Imaging Applications for Automated
Industrial Inspection and Assembly, pp. 165-170 (1979).
7. Jacobs, A. M., "First Report on the Evaluation of Gamma-
Ray Scattering Techniques for the Tomographic
Measurement of Water Two-Phase Density Distribution,"
EG&G Consultant Agreement No. K-8118-TKP-70-79, EG&G
Idaho, Idaho Falls, Idaho, July (1979).
8. Peschmann, K.R., Napel, S., Couch, J.L., Rand, R.E.,
Alei, R., Ackelsberg, S.M., Gould, R, Boyd, D.P., "High-
Speed Computed Tomography: Sytems and Performance,"
Applied Optics, Vol. 24, No. 23, pp. 4052-60 (1985).
9. Kondic, N. N., Jacobs, A. M., Anghaie, S., "Non¬
intrusive Determination of Mass Density Distribution by
Gamma-Ray Scattering: Historical Survey and
191

192
Perspectives," American Society of Mechanical Engineers-
Fluids Engineering Conference, Cavitation and Multiphase
Flow Forum, Albuquerque, New Mexico, June (1985).
10. Odeblad, E., Norhagen, A., "Electron Density in a
Localized Volume by Compton Scattering," Acta
Radiológica, Vol. 45, p. 161 (1959).
11. Lale, P. G., "The Examination of Internal Tissues Using
Gamma-Ray Scatter with a Possible Extension to Megavolt
Radiography," Physics in Medicine and Biology, Vol. 4,
p. 159 (1959).
12. Lale, P. G., "Examination of Internal Tissues by High-
Energy Scattered X-Radiation," Radiology, Vol. 90, p.510
(1968) .
13. Kondic, N. N., Hahn, O. J., "Theory and Application of
the Parallel and Diverging Radiation Beam Method in Two-
Phase Systems," Paris, France, 4th International Heat
Transfer Conference, Vol. 7, MT-1.4 (1970).
14. Farmer, F. T., Collins, M. P., "A New Approach to the
Determination of Anatomical Cross-sections of the Body
by Compton Scattering of Gamma-Rays," Physics in
Medicine and Biology, Vol. 16, No. 4, p. 557 (1971).
15. Farmer, F. T., Collins, M. P., "A Further Appraisal of
the Compton Scattering Method for Determining Anatomical
Cross-sections of the Body," Physics in Medicine and
Biology, Vol. 19, No. 6, p. 808 (1974).
16. Anghaie, S., Jacobs, A. M., Kondic, N. N., "Application
of Gamma Scattering Compton Profile for Two-Phase
Systems Densitometry," invited paper, American Nuclear
Society Winter Annual Meeting, San Francisco,
California (November, 1981).
17. Anghaie, S., "Accuracy Improvement for Gamma-Ray
Techniques in Two-Phase Measurements," Ph.D. Thesis,
Pennsylvania State University, Nuclear Engineering
(1982).
18. Kondic, N., "Density Field Determination by an External
Stationary Radiation Source Using a Kernal Techniques,"
Symposium on Polyphase Flow Measurement, ASME Winter
Annual Meeting, San Francisco (December, 1973).
19. Banerjee, S., Lahey, R. T., Jr., "Advances in Two-Phase
Flow Instrumentation," Advances in Nuclear Science and
Technology, Vol. 13, Plenum Publishing Co., New York,
New York, pp. 227-414 (1981).

193
20. Hussein, E., "A Fast Neutron Scattering Technique for
Measurement of Local Void Fraction Distribution in Two-
Phase Flow," Ph.D. Thesis, McMaster University (1983).
21. Bodette, D.E., Jacobs, A.M., "Tomographic Two-Phase Flow
Distribution Measurement Using Gamma-Ray Scattering,"
Proceedings of the Symposium on New Technology in
Nuclear Power Plant Instrumentation and Control,
Washington, D.C., November 28-30, 1984 (distributed
by Instrument Society of America).
22. Battista, J., Santón, L., Bronskill, M., "Compton-
Scatter Imaging of Transverse Sections: Corrections for
Multiple-Scatter and Attenuations," Physics in Medicine
and Biology, Vol. 22, p. 229 (1977).
23. Battista, J., Bronskill, M., "Compton-Scatter Tissue
Densitometry: Calculation of Single and Multiple-Scatter
Photon Fluences," Physics in Medicine and Biology,
Vol. 23, p. 1 (1978).
24. Kondic, N. N., Lassahn, G., "Nonintrusive Density
Distribution Measurement in Dynamic High Temperature
Systems," 24th International Instrumentation Symposium,
Albuquerque, New Mexico (May, 1978).
25. Kondic, N.N., Jacobs, A.M., Ebert, D., "Three-
Dimensional Density Field Distribution by External
Stationary Detectors and Gamma Sources Using Selective
Scattering," invited paper, 2nd International Topical
Meeting on Nuclear Reactor Thermal-Hydraulics, Santa
Barbara, California (January 11-14, 1981).
26. Tretiak, O. J., Metz, C., "The Exponential Radon
Transform," SIAM Journal of Applied Mathematics,
Vol. 39, No. 2, pp. 341-354 (1980).
27. Gullberg, G. T., Budinger, T. F., "The Use of Filtering
Methods to Compensate for Constant Attenuation in
Single-Photon Emission Computed Tomography," IEEE
Transactions on Biomedical Engineering, Vol. BME-28,
No. 2, pp. 142-157 (1981) .
28. Minerbo, G., "MENT: A Maximum Entropy Algorithm for
Reconstructing a Source from Projection Data," Computer
Graphics and Image Processing, Vol. 10, pg. 48 (1979).
29. Llacer, Jorge, "Tomographic Image Reconstruction by
Eigenvector Decomposition: Its Limitations and Areas of
Applicability," IEEE Transactions on Medical Imaging,
Vol. MI-1, No. 1, pg. 34 (July, 1982).

194
30. Herman, G. T., Image Reconstruction from Projections The
Fundamentals of Computerized Tomography, Academic Press,
New York (1980).
31. Gilbert, P., "Iterative Methods for the Three-
Dimensional Reconstruction of an Object from Projec¬
tions," Journal of Theoretical Biology, Vol. 36, pp.
105-17 (1972).
32. Joseph, P.M., "An Improved Algorithm for Reprojecting
Rays Through Pixel Images," IEEE Transactions on Medical
Imaging, Vol. MI-1, No. 3, pg. 192 (November, 1982).
33. Harms, A.A., Forrest, V.F., "Dynamic Effects in
Radiation Diagnosis of Fluctuating Voids," Nuclear
Science and Engineering, Vol. 46, pp. 408-413 (1971).
34. Hubbell, J.H., "Photon Mass Attenuation and Energy-
absorption Coefficients from 1 keV to 20 MeV,"
International Journal of Applied Radiation and Isotopes,
Vol. 33, pp. 1269-1290 (1982).
35. Anghaie, S., Finlon, K.D., Jacobs, A.M., "Optimizing and
Compton Profile Approach to Two-Phase Distribution
Measurement," Third Multiphase Flow and Heat Transfer
Symposium Workshop, Miami, Florida (April 18-20, 1983).
36. Minerbo, G.N., Sanderson, J.G., "Reconstruction of a
Source from a Few (2 or 3) Projections," LASL informal
report LA-6747-MS (March, 1977).

BIBLIOGRAPHY
Abramowitz, M., Stegun, I. A., Handbook of Mathematical
Functions, ninth edition, Dover Publications, Inc., New
York (1972).
Barth, K., Deimling, M., Fritschy, P., Mueller, E.,
Reinhardt, E.R., "Flow Measurement by Magnetic Resonance
Imaging," Application of Optical Instrumentation in
Medicine XIII - Medical Image Production, Processing,
and Display, ed. by R.H. Schneider and S.J. Dwyer,
Published by Society of Photo-Optical Instrumentation
Engineers, Vol. 535, pp. 261-267 (1985).
Gullberg, G. T., Budinger, T. F., "Three-Dimensional
Reconstruction in Nuclear Medicine Emission Imaging,"
IEEE Transactions on Nuclear Science, Vol. NS-21, p. 2
(1974).
Heinzerling, J., Schlindwein, M., "Non-Linear Techniques
in Multi-Spectral X-Ray Imaging," IEEE Transactions on
Nuclear Science, Vol. NS-27, No. 2, pg. 961 (April, 1980).
Holt, R.S., "Compton Imaging," Endeavour, New Series,
Vol. 9, No. 2, pg. 97 (1985).
Hussein, E., Banerjee, S., Meneley, D.A., "A New Fast
Neutron Scattering Technique for Local Void Fraction
Measurement in Two-Phase Flow," Thermal Hydraulics of
Nuclear Reactors, Vol. 2 (1983).
Hussein, E.M.A., "Numerical Method for Radiation Scatter
Imaging," IEEE/Seventh Annual Conference of the
Engineering in Medicine and Biology Society (1985).
King, J.D., "Coal Flowmeter Development," (CONF-840694)
From AR and TD and Suface Coal Gasification Instrumenta¬
tion Projects Contractors' Meeting, Morgantown, WV,
(June 13, 1984).
Rousseau, J.C., Czerny, J., Riegel, B., "Void Fraction
Measurement by Neutron Absorption or Scattering
Methods," OELD/NEA Specialists Meeting on Transient Two-
Phase Flow, Toronto, Canada (August, 1976).
195

BIOGRAPHICAL SKETCH
David E. Bodette was born on September 14, 1958 on
Midway Island, Hawaii. In June 1976, he graduated from Eau
Gallie High School in Melbourne, Florida. In December 1980,
he received a Bachelor of Science degree in nuclear
engineering from the University of Florida. He was awarded
an Institute of Nuclear Power Operations Fellowship in
September 1982 and in December 1982 received a Master of
Science degree in nuclear engineering from the University of
Florida.
In 1983, David Bodette was awarded a Department of
Energy Graduate Traineeship. As part of this traineeship he
spent three months at EG&G Idaho conducting experiments with
two-phase flow densitometry. He expects to receive a Doctor
of Philosophy from the University of Florida in December
1986.
196

I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality, as
a dissertation for the degree of Doctor of Philosophy.
Alan M. Jacobs
Professor of Nuclear Engineering
Sciences
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality, as
a dissertation for the degree of-Etoctorof Philosopi
Edward E. Carroll,
Professor of Nuclear
Sciences
engineering
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality, as
a dissertation for the degree of Doctor of Philosophy.
•/ ' //•• r/
â–  -v - /-â– â– '
Edward T. Dugan
Associate Professor of Nuclear
Engineering Sciences
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality, as
a dissertation for the degree of Doctor of Philosophy.
Robert J.
Professor
nrahan
Chemistry
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scppé^and quality, as
a dissertation for the degree of Doctor/óf Philosophy,
Vernon P. Roan, Jr.
Professor of Mechanical
Engineering

This dissertation was submitted to the Graduate Faculty of
the College of Engineering and to the Graduate School and
was accepted as partial fulfillment of the requirements
for the degree of Doctor of Philosophy.
December 1986
/ibhiüjf Ct'
Dean, College of Engineering
Dean, Graduate School

Internet Distribution Consent Agreement
In reference to the following dissertation:
AUTHOR: Bodette, David
TITLE: Tomographic two-phase flow measurement using Compton scattering of
gamma rays / (record number: 940973)
PUBLICATION DATE: 1986
y)A\/rd l< , as copyright holder for the aforementioned dissertation, hereby
I, _
grant specific and limited archive and distribution rights to the Board of Trustees of the University of Florida and its
agents. I authorize the University of Florida to digitize and distribute the dissertation described above for nonprofit,
educational purposes via the Internet or successive technologies.
This is a non-exclusive grant of permissions for specific off-line and on-line uses for an indefinite term. Off-line uses shall
be limited to those specifically allowed by "Fair Use" as prescribed by the terms of United States copyright legislation (cf,
Title 17, U.S. Code) as well as to the maintenance and preservation of a digital archive copy. Digitization allows the
University of Florida to generate image- and text-based versions as appropriate and to provide and enhance access using
search software.
This grant of permissions prohibits use of the digitized versions for commercialese or profit.
Signature of Copyright Holder
¿I l^oDzTr^
Printed or Typed Name of Copyright Holder/Licensee
Personal information blurred
/ f ZoQ£>
Date of Signature
Please print, sign and return to:
Cathleen Martyniak
UF Dissertation Project
Preservation Department
University of Florida Libraries
P.O. Box 117007
Gainesville, FL 32611-7007
7/1/2008