Introduction to the Mathematics of Medical Imaging: Second Edition

The front matter includes the title page, copyright page, dedication, TOC, preface, how to use this book, and Conventions

A quantitative model of a physical system is expressed in the language of mathematics. A qualitative model often precedes a quantitative model. For many years clinicians used medical x-ray images without employing a precise quantitative model. X-rays were thought of as high frequency ‘light’ with three very useful properties:

Taken together, these properties mean that using x-rays one can “see through” a human body to obtain a shadow or projection of the internal anatomy on a sheet of film [Figure 1.1(a)].

This model was adequate given the available technology. In their time, x-ray images led to a revolution in the practice of medicine because they opened the door to non-invasive examination of internal anatomy. They are still useful for locating bone fractures, dental caries, and foreign objects, but their ability to visualize soft tissues and more detailed anatomic structure is limited. There are several reasons for this. The principal difficulty is that an x-ray image is a two-dimensional representation of a three-dimensional object. In Figure 1.1(b), the opacity of the film at a point on the film plane is inversely proportional to an average of the density of the object, measured along the line joining the point to the x-ray source. This renders it impossible to deduce the spatial ordering in the missing third dimension.

Using measurements to determine the state of a system eventually reduces to the problem of solving systems of equations. In almost all practical applications, we are eventually reduced to solving systems of linear equations. This is true even for physical systems that are described by non-linear equations. As we saw in Section 1.1.1, nonlinear equations may be approximated by linear equations. A nonlinear equation is usually solved iteratively, where the iteration step involves the solution of the approximating linear system.

Linear algebra is the mathematical tool most often used in applied subjects. There are many reasons why linear equations and linear models are ubiquitous. As we will see in this chapter, there is a complete mathematical theory for systems of linear equations. Indeed, there are necessary and sufficient conditions for linear equations to have solutions and a description of the space of solutions when they exist. More pragmatically, there are efficient, generally applicable algorithms for solving systems of linear equations. This is not true even for the simplest systems of nonlinear equations. On a more conceptual level, for systems with some sort of intrinsic smoothness, a linear model often suffices to describe small deviations from a known state.

In this chapter we review some basic concepts from linear algebra, in particular the theory of systems of linear equations. This is not intended to serve as a text on finite-dimensional linear algebra but rather to situate these familiar concepts in the context of measurement problems. A more complete introduction to the theory of vector spaces is presented in Appendix A.2. As the models that arise in medical imaging are not usually finite dimensional, we briefly consider linear algebra on infinite-dimensional spaces. As a prelude to our discussion of the Fourier transform, the chapter concludes with an introduction to complex numbers and complex vector spaces.

We begin our study of medical imaging with a mathematical model of the measurement process used in x-ray tomography. The model begins with a quantitative description of the interaction of x-rays with matter called Beer's law. In this formulation the physical properties of an object are encoded in a function μ, called the attenuation coefficient. This function quantifies the tendency of an object to absorb or scatter x-rays. Using this model, idealized measurements are described as certain averages of μ. In mathematical terms, the averages of μ define a linear transformation, called the Radon transform. These measurements are akin to the determination of the shadow function of a convex region described in Section 1.2. In Chapter 6 we show that if the Radon transform of μ could be measured exactly, then μ could be completely determined from this data. In reality, the data collected are a limited part of the mathematically “necessary” data, and the measurements are subject to a variety of errors. In later chapters we refine the measurement model and reconstruction algorithm to reflect more realistic models of the physical properties of x-rays and the data that are actually collected.

This chapter begins with a model for the interaction of x-rays with matter. To explicate the model we then apply it to describe several simple physical situations. After a short excursion into the realm of x-ray physics, we give a formal definition of the Radon transform and consider some of its simpler properties. The last part of the chapter is more mathematical, covering the inversion of the Radon transform on radial functions, the Abel transform, and a discussion of Volterra integral equations.

3.1 Tomography

Literally, tomography means slice imaging. Today this term is applied to many methods used to reconstruct the internal structure of a solid object from external measurements. In x-ray tomography the literal meaning persists in that we reconstruct a three-dimensional object from its two-dimensional slices. Objects of interest in x-ray imaging are described by a real-valued function defined on ℝ³, called the attenuation coefficient. The attenuation coefficient quantifies the tendency of an object to absorb or scatter x-rays of a given energy. This function varies from point-to-point within the object and is usually taken to vanish outside it. The attenuation coefficient, like density, is nonnegative. It is useful for medical imaging because different anatomical structures have different attenuation coefficients. Bone has a much higher attenuation coefficient than soft tissue and different soft tissues have slightly different coefficients. For medical applications it is crucial that normal and cancerous tissues also have slightly different attenuation coefficients.

In this chapter we introduce the Fourier transform and review some of its basic properties. The Fourier transform is the “Swiss army knife” of mathematical analysis; it is a sturdy general-purpose tool with many useful special features. The Fourier transform makes precise the concept of decomposing a function into its harmonic components. In engineering it is used to define the power spectrum and describe many filtering operations. Because it has efficient numerical approximations, it forms the foundation for most image and signal processing algorithms.

Using its properties as a linear transformation of infinite-dimensional, normed vector spaces, we define the Fourier transform and its inverse for several different spaces of functions. Among other things we establish the Parseval formula relating the energy content of a function to that of its Fourier transform. We also study the connections between the smoothness and decay of a function and that of its Fourier transform.

In marked contrast to the Radon transform, the theory of the Fourier transform is largely independent of the dimension: The theory of the Fourier transform for functions of one variable is formally the same as the theory for functions of 2, 3, or n variables. For simplicity we begin with a discussion of the basic concepts for functions of a single variable, though in some definitions, where there is no additional difficulty, we treat the general case from the outset. The chapter ends with a recapitulation of the main results for functions of n variables.

4.1 The Complex Exponential Function

See: 2.3, A.3.1.

The building block for the Fourier transform is the complex exponential function, e^ix.

Suppose we would like to study a smooth function of one variable, but the available data are contaminated by noise. For the purposes of the present discussion, this means that the measured signal is of the form f = s + ∊n. Here ∊ is a (small) number and n is function which models the noise. An example is shown in Figure 5.1.

Noise is typically represented by a rapidly varying function that is locally of mean zero. This means that, for any x, and a large enough δ, the average

is small compared to the size of n. The more random the noise, the smaller δ can be taken. On the other hand, since s is a smooth function, the analogous average of s should be close to s(x). The moving average of f is defined to be

The Radon transform is at the center of the mathematical model for the measurements made in x-ray CT. In this chapter we prove the central slice theorem, establishing a connection between the Fourier and Radon transforms. Using it, we derive the filtered back-projection formula, which provides an exact inverse for the Radon transform. This formula is the basis of essentially all reconstruction algorithms used in x-ray CT today.

As its name suggests, the filtered back-projection formula is the composition of a one-dimensional filter with the back-projection operation, defined in Section 3.4.2. The filtration step is built out of two simpler operations: the Hilbert transform and differentiation. While differentiation is conceptually simple, it leads to difficulties when attempting to implement the filtered back-projection formula on noisy, measured data. The Hilbert transform is a nonlocal, bounded operation that is harder to understand intuitively but rather simple to implement using the Fourier representation. We analyze it in detail.

This chapter begins by reviewing the definition of the Radon transform and establishing properties of the Radon transform that are analogues of those proved for the Fourier transform in the two preceding chapters. After analyzing the exact inversion formula, we consider methods for approximately inverting the Radon transform that are relevant in medical imaging. The chapter concludes with a short introduction to the higher-dimensional Radon transform.

6.1 The Radon Transform

In Section 1.2.1 we identified ℝ × S¹ with the space of oriented lines in ℝ². The pair (t, ω) corresponds to the line

Here is the unit vector perpendicular to ω with the orientation determined by

The variable t is called the affine parameter; it is the oriented distance of the line l_t,ω to the origin.

In applications data are never collected along the whole real line or from the entire plane. Real data can only be collected from a bounded interval or planar domain. In order to use the Fourier transform to analyze and filter this type of data, we can either:

If the data are extended “by zero” then the Fourier transform is available. If the data are extended periodically then the data do not vanish at infinity, and hence their Fourier transform is not a function.

Fourier series provide a tool for the analysis of functions defined on finite intervals in ℝ or products of intervals in ℝⁿ. The goal of Fourier series is to express an “arbitrary” periodic function as a infinite sum of complex exponentials. Fourier series serve as a bridge between continuum models, like that presented in the previous chapter, and finite models, which can be implemented on a computer. The theory of Fourier series runs parallel to that of the Fourier transform presented in Chapter 4. After running through the basic properties of Fourier series in the one-dimensional case, we consider the problem of approximating a function by finite partial sums of its Fourier series. This leads to a discussion of the Gibbs phenomenon, which describes the failure of Fourier series to represent functions near to jump discontinuities. The chapter concludes with a brief introduction to Fourier series in ℝⁿ.

7.1 Fourier Series in One Dimension*

See: A.5.1, B.2.

To simplify the exposition, we begin with functions defined on the interval [0, 1]. This does not limit the generality of our analysis, for if g is a function defined on an interval [a, b], then setting f(x) = g(a+(b−a)x) gives a function defined on [0, 1] that evidently contains the same information as g.

See: A.6.1, A.5.2.

In Chapters 4 through 7, we developed the mathematical tools needed to describe functions of continuous variables and methods to analyze and reconstruct them. This chapter continues the transition from the world of pure mathematics to its application to problems in image reconstruction. In the first sections of this chapter, we imagine that our “image” is a function, f, of a single real variable, x. In a purely mathematical context, x is a real number that can assume any value along a continuum of numbers. The function also takes values in a continuum, either in ℝ, ℂ, or perhaps ℝⁿ. In practical applications we can only evaluate f at a finite set of points {x_j}. This is called sampling. As most of the processing takes place in digital computers, both the points {x_j} and the measured values {f(x_j)} are forced to lie in the preassigned, finite set of numbers known to the computer. This is called quantization. The reader is urged to review Section A.1, where these ideas are discussed in some detail.

Except for a brief discussion of quantization, this chapter is about the consequences of sampling. We examine the fundamental question: How much information about a function is contained in a finite or infinite set of samples? Central to this analysis is the Poisson summation formula. This formula is a bridge between the Fourier transform and the Fourier series. While the Fourier transform is well suited to an abstract analysis of image or signal processing, it is the Fourier series that is actually used to do the work. The reason is quite simple: The Fourier transform and its inverse require integrals over the entire real line, whereas the Fourier series is phrased in terms of infinite sums and integrals over finite intervals, both of which are eventually approximated by finite sums. Indeed, we also introduce the finite Fourier transform, an analogue of the Fourier transform for finite sequences. This chapter covers the next step from the abstract world of the infinite and the infinitesimal to the real world of the finite and discrete.

Building on the brief introduction to shift invariant filters in Section 5.1.2, this chapter discusses basic concepts in filtering theory. As noted before, a filter is the engineering term for any process that maps an input or collection of inputs to an output or collection of outputs. As inputs and outputs are generally functions of a variable (or variables), in mathematical terms, a filter is nothing more or less than a map from one space of functions to another space of functions. A large part of our discussion is devoted to linear filters, recasting our treatment of the Fourier transform in the language of linear filtering theory.

In the first part of the chapter we introduce engineering vocabulary for concepts already presented from a mathematical standpoint. In imaging applications the functions of interest usually depend on two or three spatial variables. The measurements themselves are, of necessity, also functions of time, though this dependence is often suppressed or ignored. In most of this chapter we consider filters acting on inputs that are functions of a single variable. There are three reasons for doing this: (1) The discussion is simpler; (2) this development reflects the origins of filtering theory in radio and communications; and (3) filters acting on functions of several variables are usually implemented “one variable at a time.” Most linear filters are expressed, at least formally, as integrals. The fact that higher-dimensional filters are implemented one variable at a time reflects the fact that higher-dimensional integrals are actually computed as iterated, one-dimensional integrals, that is,

In Section 9.4 we present some basic concepts of image processing and a filtering theoretic analysis of some of the hardware used in x-ray tomography.

See: A.5.2, A.6.1.

In practical applications we cannot measure an input continuously; rather it is sampled at a discrete sequence of points. The models for physical systems presented thus far are in the language of functions of continuous variables, so the question then arises as to how such models can be utilized in a real situation. In this section we analyze the transition from continuous parameter, shift invariant linear filters to their approximate realizations on finitely sampled data. For brevity we refer to this as “implementing” a shift invariant filter, with the understanding that such implementations are, of necessity, approximate. The main tool is the finite Fourier transform, and the main thing to understand is the sense in which it provides an approximation to the continuum Fourier transform. The one-dimensional case is treated in detail; higher-dimensional filters are considered briefly.

In this section x denotes an input and ℋ a shift invariant linear filter with impulse response h and transfer function ĥ. We assume that x is sampled at the discrete set of equally spaced points

The positive number τ is the sample spacing. In certain situations the physical measuring apparatus makes it difficult to collect equally spaced samples. In such situations the data are often interpolated to create an equally spaced set of samples for further analysis. This is done, even though it introduces a new source of error, because algorithms for equally spaced samples are so much simpler than those for unequal sample spacing.

In terms of the continuous parameters t and ζ, there are two different representations for ℋ: the spatial domain representation as a convolution

and the frequency domain representation as a Fourier integral

At long last we are returning to the problem of reconstructing images in x-ray tomography. Recall that if f is a function defined on ℝ² that is bounded and has bounded support, then its Radon transform, ℛ f, is a function on the space of oriented lines. The oriented lines in ℝ² are parameterized by pairs (t, ω) ∊ ℝ × S¹, with

The positive direction along l_t,ω is defined by the unit vector orthogonal to ω with det(ω ) = ±1. The Radon transform of f is the function defined on ℝ × S¹ by

Exact inversion formulæ for the Radon transform are derived in Chapter 6. These formulæ assume that ℛ f is known for all lines. In a real x-ray CT machine, ℛ f is approximately sampled at a finite set of points. The goal is to construct a function defined on a discrete set that is, to the extent possible, samples or averages of the original function f. In this chapter we see how the Radon inversion formulæ lead to methods for approximately reconstructing f from realistic measurements. We call such a method a reconstruction algorithm.

In this chapter we derive the reconstruction algorithms used in most modern CT scanners. After a short review of the basic setup in x-ray tomography, we describe the principal designs used in modern CT scanners. This, in turn, determines how the Radon transform of the attenuation coefficient is sampled. We next derive algorithmic implementations of the filtered back-projection formula for parallel beam, fan beam, and spiral CT scanners. A large part of our discussion is devoted to an analysis of the filtration step and its efficient implementation.

In the previous chapter we derived finite algorithms to approximately reconstruct a function of two variables with bounded support from finitely many samples of its Radon transform. For a variety of reasons this is still a highly idealized situation. In this chapter we analyze these algorithms with a more realistic model for the measurement process. The first issue we address is the fact that the x-ray beam has a finite (two-dimensional) width. There is a simple linear model for this effect as weighted averaging in the affine parameter. More careful investigation reveals that this is a nonlinear phenomenon that leads to the nonlinear partial volume effect. Averaging in the affine parameter is a form of lowpass filtering, which is, in turn, important for the analysis of the aliasing that results from sampling the Radon transform.

We next derive the total point spread function (PSF) for the measurement and reconstruction process first without and then with sampling. Once the measurements are sampled in the affine parameter, the reconstruction process is no longer shift invariant, so it is not described by a single point spread function. The width of the central peak of the PSF provides an indication of the resolution available in the reconstructed image. We analyze the effects of the beam width and sample spacing on the shape of the PSF. Using the PSF, we consider the consequences for the reconstructed image of various sorts of measurement errors. At the end of the chapter we consider beam hardening; this results from the fact that the x-ray beam is not monochromatic. Beam hardening is another fundamentally nonlinear phenomenon.

12.1 The Effect of a Finite Width X-Ray Beam

Up to now, we have assumed that an x-ray beam is just a line with no width and that the measurements are integrals over such lines. What is really measured is better approximated by averages of such integrals. We now consider how the finite width of the x-ray beam affects the measured data. Our treatment closely follows the discussion in [113].

Algebraic reconstruction techniques (ARTs) are techniques for reconstructing images that have no direct connection to the Radon inversion formula. Instead these methods use a linear model for the measurement process so that the reconstruction problem can be posed as a system of linear equations. Indeed the underlying mathematical concepts of ART can be applied to approximately solve many types of large systems of linear equations. This chapter contains a very brief introduction to these ideas. We present this material for historical reasons and to expose the reader to other approaches to medical image reconstruction. An extensive discussion of these methods can be found in [52]. The material presented in this chapter in not used in the sequel.

As before, boldface letters are used to denote vector or matrix quantities, while the scalar entries of a vector or matrix are denoted in regular type with subscripts. For example, r is a matrix, {r_i} its rows (which are vectors), and r_ij its entries (which are scalars).

13.1 Algebraic Reconstruction

The main features of the Radon transform of interest in ART are as follows: (1) The map f → ℛ f is linear, and (2) for a function defining a simple object with bounded support, ℛ f has a geometric interpretation. The first step in an algebraic reconstruction technique is the choice of a finite collection of basis functions,

Certain types of a priori knowledge about the expected data and the measurement process itself can be “encoded” in the choice of the basis functions. At a minimum, it is assumed that the attenuation coefficients we are likely to encounter are well approximated by functions in the linear span of the basis functions.

14.1 Introduction

Nuclear magnetic resonance (NMR) is a subtle quantum mechanical phenomenon that has played a major role in the revolution in medical imaging over the last 30 years. Before being used in imaging, NMR was employed by chemists to do spectroscopy, and remains a very important technique for determining the structure of complex chemical compounds like proteins. There are many points of contact between these technologies, and problems of interest to mathematicians. Spectroscopy is an applied form of spectral theory. NMR imaging is connected to Fourier analysis, and more general Fourier Integral Operators. The problem of selective excitation in NMR is easily translated into a classical inverse scattering problem, and in this formulation, is easily solved.

In this chapter, we give a rapid introduction to nuclear magnetic resonance imaging. While the physical basis for this technique is purely quantum mechanical, i.e., the fact that the protons in hydrogen atoms are spin-1/2 particles, there is an entirely classical model that is adequate to understand most applications of NMR in medical imaging. This is called the Bloch phenomenological equation, and it forms the foundation of our presentation. X-ray tomography is based upon a much simpler physical phenomenon: the scattering and absorption of x-ray photons, and therefore offers very few mechanisms for creating contrast between different tissue types, and essentially no mechanism for imaging metabolic processes. NMR is based on subtle and complex physical phenomena, which allows for a great variety of different contrast mechanisms.

In our treatment, little attention is paid to either NMR spectroscopy, or the quantum description of NMR. Those seeking a more complete introduction to these subjects should consult the monographs of Abragam, [1], Levitt [84], or Ernst, Bodenhausen and Wokaun, [38], for spectroscopy, and those of Callaghan, [17], Haacke, et al. [50], and Bernstein, et al. [9], for imaging. All six books consider the quantum mechanical description of the these phenomena.

Up to this point, we have only considered deterministic models. These are models where a known input produces a definite output. We now begin to discuss probability theory, which is the language of noise analysis. But what is noise? There are two essentially different sources of noise in the mathematical description of a physical system or measurement process. The first step in building a mathematical model is to isolate a physical system from the world in which it sits. Once such a separation is fixed, the effects of the outside world on the state of the system are often modeled, a posteriori, as noise. In our model for CT imaging it is assumed that every x-ray that is detected is produced by our x-ray source. In reality, there are many other sources of x-rays that might also impinge on our detectors. Practically speaking, it is not possible to give a complete description of all such external sources, but sometimes it is possible to describe them probabilistically.

Many physical processes are inherently probabilistic, and that is the second source of noise. In CT imaging the interaction of an x-ray “beam” with an object is a probabilistic phenomenon. The beam is, in fact, a collection of discrete photons. Whether or not a given photon, entering an object on one side, reemerges on the other side, traveling in the same direction, depends on the very complicated interactions this photon has with the microscopic structure of the object it passes through. These individual interactions cannot be modeled in a useful way. If μ is the attenuation coefficient and I (t) is the incident flux of photons at t, then Beer's law says that the change in the flux, I (t + Δt) − I (t), over a small distance Δt is

or

This can be interpreted as the statement that an incident photon has probability (1 − μ Δt) of being emitted. Beer's law describes the average behavior of a photon and is useful only if the x-ray beam is composed of a large number of photons.

In this chapter we consider several applications of the theory of probability to medical imaging. We begin with probabilistic models for x-ray generation and detection and a derivation of Beer's law. In the next section we analyze the propagation of noise through the filtered back-projection algorithm. This leads to relations between the signal-to-noise ratio and the sample spacing. Using the connection between the sample spacing and spatial resolution, this, in turn, gives a relationship between the dosage of radiation and the resolution at a fixed signal-to-noise ratio. In the last section we briefly describe positron emission tomography and the maximum likelihood algorithm for reconstructing images. This algorithm is a probabilistic alternative to filtered back-projection that may give superior results if the measurements are very noisy.

16.1 Applications of Probability Theory to X-Ray Imaging

In this section we present simple probabilistic models for an x-ray detector and a source-detector pair. We also give a probabilistic derivation of Beer's law.

16.1.1 Modeling a Source-Detector Pair

Suppose that we have a Poisson source with intensity and a Bernoulli detector. This means that each incident photon has probability p, 0 ≤ p ≤ 1 of being detected and each detection event is independent of any other. What is the statistical description of the output of such a source-detector pair? The probability of observing k photons, given that N photons arrive, is a conditional probability:

To model noise in measurement and filtering processes requires concepts more general than that of a random variable. This is because we need to discuss the results of passing a noisy signal through a linear filter. As was the case in chapter 9, it is easier to present this material in terms of functions of a single variable, though the theory is easily adapted to functions of several variables. After giving the basic definitions, we consider some important examples of random processes and analyze the effect of a linear shift invariant filter on random inputs. The chapter concludes with an analysis of noise in the continuum model for filtered back-projection. Our discussion of random processes is very brief, aimed squarely at the goal of analyzing the effects of noise in the filtered back-projection algorithm.

17.1 Random Processes in Measurements

To motivate this discussion, we think of the familiar example of a radio signal. A radio station broadcasts a signal s_b(t) as an electromagnetic wave. A receiver detects a signal s_r(t), which we model as a sum

Here F is a filtering operation that describes the propagation of the broadcast signal. For the purposes of this discussion we model F as attenuation; that is, F(s_b) = λs_b for a 0 < λ < 1. The other term is noise. In part, the noise records “random” aspects of the life history of the broadcast signal that are not modeled by F. Beyond that it is an aggregation of other signals that happen to be present at the carrier frequency of the broadcast signal. The existence of the second part is easily verified by tuning an (old) radio to a frequency for which there is no local station. Practically speaking, it is not possible to give a formula for the noise. Because we cannot give an exact description of the noise, we instead describe it in terms of its statistical properties.

We begin with the assumption that the noise is a bounded function of t. What can be done to specify the noise further? Recall the ensemble average definition of the probability of an event as the average of the results of many “identical experiments.” In the radio example, imagine having many different radios, labeled by a set . For each α ∊ we let s_r,a(t) be the signal received by radio α at time t. The collection of radios is the sample space. The value of s_r,a at a time t is a function on the sample space—in other words, a random variable. From the form given previously we see that

In applied subjects, mathematics needs to be appreciated in three rather distinct ways: (1) in the abstract context of perfect and complete knowledge generally employed in mathematics itself; (2) in a less abstract context of fully specified, but incompletely known functions—this is the world of mathematical approximation; and (3) in a realistic context of partially known functions and noisy, approximate data, which is closer to the real world of measurements. With these different perspectives in mind, we discuss some of the mathematical concepts underlying image reconstruction and signal processing. The bulk of this material is usually presented in undergraduate courses in linear algebra, analysis, and functional analysis. Instead of a giving the usual development, which emphasizes mathematical rigor and proof techniques, we present this material from an engineering perspective. Many of the results are proved in exercises, and examples are given to illustrate general phenomena. This material is intended to provide background material and recast familiar material in a more applied framework; it should be referred to as needed.

A.1 Numbers

We begin by discussing numbers, beginning with the abstract concepts of numbers and their arithmetic properties. Representations of numbers are then considered, leading to a comparison between abstract numbers and the way numbers are actually used in computation.

A.1.1 Integers

Mathematicians think of numbers as a set that has two operations, addition and multiplication, that satisfy certain properties. The mathematical discussion of this subject always begins with the integers. We denote the set of integers by ℤ and the set of positive integers (the whole or natural numbers) by ℕ. There are two operations defined on the integers: addition, +, and multiplication, ×. Associated to each of these operations is a special number: For addition that number is 0 it is defined by the property

This appendix contains some of the basic definitions and results from analysis that are used in this book. Many good treatments of this material are available; for example, [27], [111], or [121].

B.1 Sequences

A sequence is an ordered list of objects. As such it can be described by a function defined on a subset of the integers of the form {M, M + 1, …, N − 1, N}. This set is called the index set. If both M and N are finite, then the sequence is finite. If M = −∞ and N is finite or M is finite and N = ∞, then it is an infinite sequence. If the index set equals ℤ, then the sequence is bi-infinite. In this section and throughout most of the book, the term sequence usually refers to an infinite sequence. In this case the index set is usually taken to be the positive integers ℕ. For example, a sequence of real numbers is specified by a function

With this notation the nth term of the sequence would be denoted x(n). It is not customary to use functional notation but rather to use subscripts to label the terms of a sequence. The nth term is denoted by x_n and the totality of the sequence by < x_n >. This distinguishes a sequence from the unordered set consisting of its elements, which is denoted {x_n}.

Given a sequence < x_n >, a subsequence is defined by selecting a subset of {x_n} and keeping them in the same order as they appear in < x_n >. In practice, this amounts to defining a monotone increasing function from ℕ to itself. We denote the value of this function at j by n_j. In order to be monotone increasing, n_j < n_j+1, for every j in ℕ. The jth term of the subsequence is denoted by x_nj, and the totality of the subsequence by < x_nj >. As an example, consider the sequence x_n = (−1)ⁿn; setting n_j = 2 j defines the subsequence x_nj = (−1)^2j2j.

Definition B.1.1. A sequence of real numbers, < x_n >, has a limit if there is a number L such that, given any ∊ > 0, there exists an integer N > 0, such that

The back matter includes bibliography and index