A surprising radon transform result and its ... - Semantic Scholar
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
1039
A Surprising Radon Transform Result and Its Application to Motion Detection Thomas L. Marzetta, Senior Member, IEEE, and Larry A. Shepp
Abstract—An elliptical region of the plane supports a positivevalued function whose Radon transform depends only on the slope of the integrating line. Any two parallel lines that intersect the ellipse generate equal line integrals of the function. We prove that this peculiar property is unique to the ellipse; no other convex, compact region of the plane supports a nonzero-valued function whose Radon transform depends only on slope. We motivate this problem by considering the detection of a constant-velocity moving object in a sequence of images, in the presence of additive, white, Gaussian noise. The intensity distribution of the object is known, but the velocity is only assumed to lie in some known set, for example, an ellipse or a rectangle. The object is to find a space-time linear filter, operating on the image sequence, whose minimum output signalto-noise ratio (SNR) for any velocity in the set is maximized. For an ellipse (and its special cases, the disk and the linesegment) the special Radon transform property of the ellipse enables us to obtain a closed-form, analytical solution for the minimax filter, which significantly outperforms the conventional three-dimensional (3-D) matched filter. This analytical solution also suggests a constrained minimax filter for other velocity sets, obtainable in closed form, whose SNR can be very close to the minimax SNR. Index Terms—Convex set, motion detection, Radon transform.
I. INTRODUCTION
A
CHALLENGING problem in image sequence analysis is to detect a moving object in the presence of additive noise [1], [3], [4], [14]–[18]. It is convenient and often reasonable to make several assumptions: that the object is moving with constant velocity, that its image intensity distribution does not change with time, and that the noise is zero mean Gaussian that is uncorrelated over both space and time. We denote the where represents time noisy image sequence by is Cartesian position within the image frame and (1) is the intensity distribution of the object, is the velocity of the object expressed in some consistent set of units, for example meters/s or pixels/frame, is the position of the object at reference time , and is additive white noise. By assumption, the noise variance is so large that the object cannot be detected
where
Manuscript received June 24, 1997; revised April 3, 1998. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Steven D. Blostein. T. L. Marzetta is with Bell Laboratories, Lucent Technologies, Murray Hill, NJ 07974-0636 USA (e-mail: [email protected]). L. A. Shepp is with the Department of Statistics, Rutgers University, New Brunswick, NJ 08903 USA (e-mail: [email protected]). Publisher Item Identifier S 1057-7149(99)06116-3.
reliably merely by thresholding the image sequence, or by applying a two-dimensional (2-D) (spatial) matched filter to each image frame. In general we need a more sophisticated three-dimensional (3-D) detection scheme that operates over multiple pixels and multiple frames of the image sequence. A. Four Motion Detection Problems The nature of the detection problem depends on the extent of prior information that is available concerning the object of interest. In almost any practical problem the reference position , is unknown. The velocity may or may of the object, not be known, and the object intensity distribution may or may not be known. It is instructive to consider briefly the four cases, each of which requires a distinct detection strategy. The fourth case, an unknown velocity and a known object intensity distribution, is the most interesting case, as well as the focus of this paper. 1) Velocity Known; Object Intensity Distribution Known: Here we know everything about the signal except for the reference position. Correlating the noisy image sequence with a replica of the signal to be detected for all possible generates a set of sufficient statistics. The actual detection operation compares the largest value of the correlation with a threshold; this constitutes a generalized likelihood ratio test [20] for the presence or absence of the moving object. It is necessary to limit the temporal duration of the correlation operation if only because the object remains in the field of view for a finite length of time. The 3-D (i.e., space-time) correlation operation is equivalent to applying the so-called 3-D matched filter [16], whose impulse response is (2) is a finite duration temporal window, for example where a rectangle or a raised cosine. When applied to the noisy image sequence the filter produces an image sequence containing a moving object having the original velocity, and with an improved signal-to-noise ratio (SNR). The 3-D matched filter operates in two ways: it implements a spatial matched filter, and it integrates the image sequence along lines in the space that coincide with the trajectory of the moving object. The signal-to-noise improvement is proportional to the duration of the temporal window. 2) Velocity Known; Object Intensity Distribution Unknown: is unknown, and that the image sequence Suppose that is observed over the time interval, [0, ]. It is straightforward to show that, with knowledge only of the velocity, a sufficient
1057–7149/99$10.00 1999 IEEE
1040
Fig. 1. Bank of 3-D matched filters can cover a region of possible object velocities. Filters achieve their maximum SNR at the grid points; velocity mismatch occurs between grid points.
statistic is generated by integrating the image sequence along the trajectory of the moving object. This sufficient statistic is a single image that contains the object with an enhanced SNR that is proportional to
Therefore the original detection problem reduces to detecting an unknown object in a single noisy image. To proceed any . further it is necessary to assume something about For example if the extent of the object were known, one with a disk or rectangle whose size could correlate , is comparable to that of the object. The generation of and the spatial correlation operation together are equivalent to is replaced applying the 3-D matched filter (2), where by a disk or rectangle within the formula. 3) Velocity Unknown; Object Intensity Distribution Unknown: Situations where the velocity is unknown occur more frequently than where the velocity is known. First consider and the velocity are unknown. the case where both Conceptually one could repeatedly apply the 3-D matched filter for all possible velocities—an impractical scheme. In practice a bank of 3-D matched filters has been proposed as in Fig. 1, whose design velocities span the possible range of velocity of the object [1], [15], [17]. Given the severe computational burden of implementing a 3-D filter it is desirable to maximize the efficiency of the filterbank to achieve reliable detection with as few filters as possible by optimizing the temporal window and the design velocities. When the object velocity falls between velocity-grid points in the filterbank velocity mismatch occurs, and the 3-D matched filter is no longer optimal; in fact, increasing the temporal duration of the filter indefinitely eventually results in a reduced SNR, because then the filter integrates mostly noise. The optimization of a bank of 3-D matched filters can be difficult because there is no general closed-form expression for SNR improvement when there is velocity mismatch. An alternative approach handles velocity uncertainty by designing the filter to have an adequate SNR over a set of velocities—for example, a square or a hexagon—rather than
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Fig. 2. Filterbank where each filter yields the stipulated SNR uniformly for every velocity inside of its tile.
to have a high SNR at merely one velocity [11]–[13]. One can design a bank of such filters for a larger velocity set by tiling the region of interest in the velocity plane as in Fig. 2, such that each tile corresponds to one of the filters. For the case where the velocity is unknown and the object intensity distribution is unknown, one can require that a 3-D filter have minimum white-noise sensitivity subject to the constraint that it pass all moving objects, whose velocity is (that corresponds to one of the tiles contained in a set in Fig. 2) with unity gain [11], [12]. This design procedure yields a constant SNR over the velocity set. The resulting filter is an example of a fan filter—originally used in geophysical applications to filter wavefields on the basis of plane wave slowness [6], [7], [9]. Since the fan filter has unity gain over the 3-D frequency support of the moving object, it does not destroy any information. For elementary velocity sets there are closed-form expressions for the frequency response of the filter and for its SNR improvement. 4) Velocity Unknown; Object Intensity Distribution Known: In certain cases, the object intensity distribution is known, but the velocity of the object is unknown. This can occur, for example, when the angular extent of the object is much smaller than the angular resolution of the sensor, in which case the effective object intensity distribution is equal to the known point spread function of the sensor. An intuitively appealing approach to this problem would cascade a spatial matched filter (e.g., matched to the known object intensity distribution) with the 3-D fan filter mentioned in Section I-A3. It turns out, however, that this combination can be improved upon. Instead, without imposing any a priori structure one can formulate a minimax design problem for the filter where, for all velocities in the set , the minimum SNR is maximized [13]. A general analytical solution to the minimax problem was not found, but for the case where is an elliptical region a closed-form analytical solution was found whose existence depends on the set supporting a peculiar function. B. Peculiar Radon Transform Property In the case of an elliptical velocity set, the key to the analytical solution for the minimax filter is the existence of a positive Lagrange multiplier function on the elliptical region whose Radon transform depends only on the slope of the
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
1041
sequence yields the best SNR in a minimax sense for all velocities in the set, . The noise variance at the output of the . The result of filter is equal to filtering the moving object is again a moving object that has the same velocity, but a different object intensity distribution. We want to maximize the amplitude of the filtered moving object, which is proportional to the filter gain, at the arbitrary point, . (Shifting this point would merely shift the optimized filter impulse response.) The SNR at the output of the filter is a function of the unknown velocity, SNR Fig. 3. Convex elliptical region supports a positive-valued function such that all parallel lines that intersect the ellipse generate equal line integrals; the Radon transform depends only on the slope of the integrating line.
integrating line; any two parallel lines that intersect the ellipse generate equal line integrals of the function (see Fig. 3). It is somewhat surprising that such a function exists. It is even more surprising, as we show in this paper, that the elliptical region is the only convex compact set that supports such a function. C. Organization of Paper Section II formulates the design of the minimax filter, which is equivalent to an infinite dimensional quadratic programming problem. We apply the Kuhn–Tucker conditions to obtain the optimal filter in terms of a nonnegative-valued Lagrange , defined on the velocity set multiplier function, which, in turn, satisfies a 2-D linear integral equation. While a general solution to the integral equation eludes us, we show that an analytical solution can be obtained if a function exists that has the peculiar Radon transform property mentioned in Section I-B. There is such a function when is an ellipse. The minimax filter for the ellipse suggests a constrained minimax filter for other sets, such as the rectangle or the hexagon that, although suboptimal, can be obtained analytically, and whose performance approaches that of the optimal minimax filter. Section III compares the performance of the minimax filter, the constrained minimax filter, and the optimized 3-D matched filter for some cases of interest. Section IV provides a constructive proof that the ellipse is the only convex, compact set that supports a nonzero-valued function whose Radon transform depends only on the slope of the integrating line. II. MINIMAX 3-D FILTER
(3) where the numerator is the square of the filter gain as a function of velocity. The optimal filter is the solution to the following minimax optimization problem: SNR This minimax problem is equivalent to an infinite dimensional quadratic programming problem where we minimize the noise variance at the output of the filter subject to a linear inequality constraint on the filter gain, for all possible velocities, subject to
(4) where is an arbitrary constant (changing merely scales the amplitude of the filter without changing the SNR). It is convenient to reformulate the minimax problem in the 3-D frequency domain. The 3-D frequency response of the filter is
where we use the geophysical convention with respect to the signs in the exponent. The 3-D Fourier transform of the moving object is
A. Problem Statement We observe the noisy image sequence, (1), where the object , is known, the reference position intensity distribution, , is unknown, and the velocity of the of the object, , is unknown. However the velocity is assumed object, to lie in some known convex compact set, . The noise, , is zero mean, Gaussian, white with known spectral . density, We seek a 3-D linear shift-invariant filter having impulse , which when applied to the noisy image response,
where is the Dirac delta function, and is the 2-D Fourier transform of the object intensity distribution
1042
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Thus, the moving object occupies a plane in the 3-D frequency space that passes through the origin
(7)
(5) With the above notation the minimax problem (4) becomes (8) subject to
(9)
(6) A feasible solution (e.g., a finite energy filter that satisfies the constraints) always exists: an example is the cascade of the spatial matched filter and the fan filter mentioned in Section I-A4. For mathematical convenience, we have formulated the detection problem under the assumption that the noisy image “sequence” is measured continuously with respect to space and time, and the 3-D filter is a continuous parameter filter rather than a sampled data filter. In the case of a diffraction-limited optical sensor, the modulation transfer function is spatially bandlimited, and so is the object intensity distribution [2]; it follows from the relation (5) that for all finite velocities, the moving object is temporally bandlimited as well. Under these conditions a spatially and temporally lowpass filtered version of the continuous parameter noisy image sequence could be sampled at the Nyquist rate without loss of information, and the discrete filter would merely be a sampled version of the continuous filter. B. Kuhn–Tucker Conditions The inequality constraints involved in the optimization (6) can be converted into equality constraints by subtracting a nonnegative function of velocity from the left-hand side of the constraint; in turn this nonnegative function is equal to the . We include square of an unconstrained function, the resulting equality constraints in the minimization through the method of Lagrange multipliers. Since there are an infinite number of constraints we use a Lagrange multiplier function, . The resulting minimization problem with respect , , and becomes to
The calculus of variations gives three conditions:
The first condition (7) gives the optimal filter in terms of the Lagrange multiplier function, and the second condition (8) is merely the equality constraint for every velocity. The substitution of (7) into (8), and the integration over , yields a 2-D integral equation for the Lagrange multiplier function
(10) appears in the integral equation, Note that the function and are also mutually dependent through the and that third condition (9). The quantitative interpretation of the Lagrange multiplier is that it is the derivative of the optimized cost (e.g., the noise variance at the output of the filter) with respect to the righthand side of the constraint, [19]. As increases, implying an increase in filter gain, the noise variance at the output of the optimized filter cannot decrease. Consequently, the Lagrange multiplier function must be nonnegative. We combine this fact with (9) to obtain the following Kuhn–Tucker conditions [10]:
if
then
if
then
(11) .
is equal to zero, the filter gain is For velocities where exactly equal to and the filter achieves the minimax SNR is not equal to zero, exactly, while for velocities where the filter gain exceeds and the SNR exceeds the minimax SNR. The implications of the Kuhn–Tucker conditions are the following. must be 1) The solution of the integral equation for nonnegative. is positive the filter yields 2) For all velocities where exactly its minimax SNR. 3) For all velocities where the SNR exceeds the minimax must be equal to zero. SNR, In short, constraints are active where the Lagrange multiplier function is greater than zero. We do not have an analytical solution to the integral equation in combination with the Kuhn–Tucker conditions for the general case. However, in the next two subsections we derive some general properties of the optimal filter, and we actually obtain a closed-form analytical solution for the case is an ellipse. where
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
1043
Then one can verify by direct substitution that the solutions in the case of the translated set are
Fig. 4. Support function for a convex, compact set; distance between the tangent line and the origin.
f ( )
is the signed
C. General Properties of Minimax Filter The minimax filter (7) is equal to the product of the spatial , and a 3-D filter that consists of a linear matched filter, superposition of impulsive planes that are defined by (5). As a result, the spatial frequency support of the minimax filter is identical to that of the object intensity distribution, and , the union of the impulsive planes for a particular determines the temporal frequency support. We denote the , temporal frequency support by for some In turn this set can be conveniently expressed in terms of the support function [8] for the convex velocity set, (12) As shown in Fig. 4, the support function is equal to the signed distance between a tangent line and the origin. It follows easily that (13) Thus the convex set is defined by a pair of parallel tangent lines for every angle . We express spatial frequency in polar coordinates, (14) It follows that the temporal frequency support for the minimax filter is the continuous interval
The shift of the Lagrange multiplier function changes the frequency response (7) of the filter nontrivially, however a direct calculation of the energy of the filter shows that its noise sensitivity, and therefore the minimax SNR, is unchanged. If one is designing a bank of minimax filters using identical tiles, then only one optimization problem needs to be solved because of this property, and all of the filters yield the same minimax SNR. Property 2—For a Fixed Shape and Orientation of the Velocity Set, the Minimax SNR Is Inversely Proportional to the and Square-Root of the Area of the Set: Suppose that satisfy the integral equation and the Kuhn–Tucker conditions for some velocity set, . Let be the velocity set obtained by expanding (in both dimensions) the set by the factor , which increases the area by the factor ,
It can be shown that the solutions in the case of the expanded are set
A calculation of the energy of the resulting minimax filter discloses that the minimax SNR (in units of power/power) is decreased by the factor . For example, to increase the minimax SNR by 3 dB (i.e., doubling the SNR), we have to decrease the linear size of the velocity set by a factor of two. As a result the velocity-area of each tile is reduced by a factor of four, or equivalently the number of filters has to increase by a factor of four to cover the same region of interest. D. Closed-Form Solution When
Is Elliptical
In general we can replace the integral equation involving by an equivalent integral equation involving the Radon . The inner integral of (10) is related to the transform of Radon transform of the Lagrange multiplier function; with the polar representation of spatial frequency, (14) we have
(15) Without actually solving for the minimax filter, we can deduce two important properties that greatly simplify the task of designing a filterbank. Property 1—The Minimax SNR Is Invariant to Translations of the Velocity Set: Suppose that for the set , we have and that satisfy the integral equation (10) solutions for and the Kuhn–Tucker conditions (11). Let be the velocity set , i.e., obtained by translating the set by the amount
where
is the Radon transform of
1044
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
After substituting this expression into (10) the integral equation becomes
The calculated SNR is consistent with Property 2 of Section IIC. We show in Section IV that the ellipse is the only set that has the peculiar Radon transform property that we exploited here to obtain an analytical solution for the minimax filter.
(16) Note that the velocity affects the integral only through the argument of the Radon transform. If we can find a nonnegative Lagrange multiplier function whose Radon transform has no -dependence, for all values of such that the integrating line intersects the velocity set, then a scaled version of automatically satisfies (16) with . For example suppose that is the elliptical region, . Then it can be verified by elementary calculation that the following positive function:
E. Constrained Minimax Filter for Arbitrary Convex Velocity Sets
(17) has the Radon transform
For cases where the velocity set is not an ellipse, one could always obtain the minimax filter by numerical means. To that end, a dual optimization problem is more tractable than the primal optimization problem (6). It can be shown that the dual problem involves choosing the nonnegative-valued Lagrange multiplier function to maximize a convex quadratic functional, with no other constraints. However we do not pursue the numerical approach in this paper. Instead, we obtain a useful suboptimal filter analytically. For the case where is an ellipse, note that, within its passband, the frequency response of the minimax filter (20) is independent of . This motivates us to consider, in general, a filter of the form
otherwise
otherwise (18) where
is the support function, (12), (19)
depends only on spatial frequency, and where where again spatial frequency is expressed in polar coordinates. The SNR of this filter is the same for all velocities within . The that maximizes the SNR gives (to within a choice of multiplicative constant)
After substituting (17) into (7) and using (18), we obtain the closed form expression for the minimax filter (to within an arbitrary multiplicative constant)
(22) otherwise. We refer to this as the constrained minimax filter. The resulting SNR is
otherwise. (20)
SNR
The frequency support of the filter is an elliptical cone. The is no cause for alarm. In the singularity at , vicinity of the origin the factor, , is both absolutely- and square-integrable with . The presence of this singularity is inturespect to itively reasonable: the frequency support plane (5) associated with every velocity in passes through the origin, so it makes sense to concentrate the filter energy near the origin rather than to spread the filter energy evenly throughout the filter support. The resulting SNR is the same for all velocities inside of the ellipse:
(23) It is easily shown that the two properties of Section II-C apply to the constrained minimax filter as well. We later show that the performance of this filter can be very close to that of the optimal minimax filter. III. PERFORMANCE RESULTS In what follows we assume that the object intensity distribution is Gaussian-shaped:
SNR (21)
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
1045
(24) is the integrated intensity of the where object, and represents the size of the object. This is a common model for the point-spread function of an optical sensor. Additionally it turns out that this object intensity distribution allows us to obtain performance results analytically. A. Disk-Shaped Velocity Set: Minimax Filter Versus 3-D Matched Filter Here we assume a disk-shaped velocity set of radius centered at the origin, and we compare the performance of the minimax filter with that of an optimized 3-D matched filter. 1) Minimax Filter: The minimax filter is given by (20), , and the minimax SNR is given evaluated for by (21), easily evaluated analytically for the Gaussian object intensity distribution (24). We can also evaluate the SNR for velocities outside of the disk, which is less than the minimax SNR, SNR
. (25) 2) 3-D Matched Filter: The 3-D matched filter (2) is characterized by its design velocity and by its temporal window. Its 3-D frequency response, when tuned for zero-velocity (the center of ), is
Because of the circular symmetry of , the resulting SNR depends only on the magnitude of the velocity, and it decreases monotonically with the speed. We optimize the filter in a to maximize the SNR on minimax sense by choosing the boundary of the velocity disk [specifically at the point, ]. The result of the optimization is
Fig. 5. Filterbank based on hexagonal tiling. Each tile is covered by a constrained minimax filter that provides uniformly good SNR. The SNR of the filters is within 0.4 dB of the minimax SNR.
improvement factor is (0.5 dB).1 Equivalently, for the same guaranteed SNR, the minimax filter can cover a velocity disk whose radius is 2.2 times greater than that of the 3-D matched filter. The improved performance of the minimax filter relative to the 3-D matched filter is even more dramatic for an ellipse. The 3-D matched filter has to be designed for the major axis , and its minimum SNR within the of the ellipse, say ellipse is inversely proportional to . In contrast, the minimax remains SNR increases without limit as decreases while constant. B. Designing a Filterbank: Constrained Minimax Filters versus 3-D Matched Filters Suppose that the possible range of velocity of the object is so great that the required SNR cannot be realized by a single filter, necessitating the design of a bank of filters. We compare two approaches. The first approach tiles the velocity plane with hexagons, with each tile covered by the constrained minimax filter (22). The second approach uses a bank of optimized 3-D matched filters whose design velocities are on a hexagonal grid. 1) Constrained Minimax Filters for Hexagonal Tiles: Fig. 5 shows a hexagonal tiling of the velocity plane, with the halfspacing between the centers of adjacent tiles equal to . The support function for a hexagon centered at the origin is mod . An evaluation of (23) with this support function gives the resulting SNR for any of the constrained minimax filters when the velocity is in the corresponding tile; this expression is then solved for the required tile size:
The SNR resulting from the optimized 3-D matched filter is SNR
SNR or
(26) 3) Comparison: Comparing the SNR’s of the two filters, (25) and (26), we find that the minimax filter performs better than the 3-D matched filter everywhere in the velocity disk: on the boundary of the disk the SNR improvement factor is (3.5 dB), while at the center of the disk the
SNR
(27)
It would be interesting to compare the SNR for the constrained minimax filter with the SNR for the ideal minimax 1 An even greater improvement factor was given in [13], however a different object intensity distribution than (24) was assumed, and the 3-D matched filter was not truly optimized.
1046
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Fig. 6. Performance of minimax filter for hexagon is bounded by that of minimax filters for outer disk and inner disk.
filter, but this would require a numerical optimization. Instead we can obtain upper and lower bounds on the minimax SNR, equal to the minimax SNR’s for the two disks shown in Fig. 6, , and the inner disk the outer disk having a radius of having a radius of : SNR A comparison of the upper bound with (27) shows that the optimal minimax filter for the hexagon would have an SNR (0.42 dB) that is at most a factor of greater than that of the constrained minimax filter. The ideal minimax filter could reduce the size of the filterbank by, at most, 20%. 2) 3-D Matched Filters on a Hexagonal Velocity Grid: Fig. 7 illustrates a bank of 3-D matched filters arranged on a hexagonal velocity grid. The half-spacing between adjacent grid points is . Each filter has a disk-shaped region of , in which the filter achieves coverage, of radius the stipulated SNR. Of necessity, the disks have to be overlapped in the manner shown. Equation (26) relates to the minimum SNR that is achieved by the filterbank; this expression, in turn, gives the minimum SNR in terms of , SNR
Fig. 7. Bank of 3-D matched filters is less efficient than bank of constrained minimax filters shown in Fig. 5. Each filter achieves the stipulated SNR on the boundary of the disk, with maximum SNR at the center of the disk. Disks must overlap to compensate for velocity mismatch.
parallel lines that intersect the set generate equal line integrals of the function. We now prove that this property is unique to the elliptical region. A. Desired Radon Transform Property be a convex, compact set in the plane, and let be a function defined on . The Radon transform of the function is the line integral of the function, for the , integrating line, Let
(29) The Radon transform vanishes for outside of the support in, where is the support function, terval, (12). , whose The object is to find a nonzero-valued function, for all values of Radon transform is independent of within the support interval. If such a function exists, then it is straightforward to infer its Radon transform. To do so, we integrate both sides of (29) with respect to , as follows:
or SNR
(28)
3) Comparison: Upon comparing (27) with (28) we find that the constrained minimax filterbank is considerably more efficient than the 3-D matched filterbank. To achieve the same minimum SNR with 3-D matched filters we must have
Using the 3-D matched filter rather than the constrained times as many minimax filter would require filters. IV. NEW RADON TRANSFORM RESULT In the previous section, we found a closed-form analytical solution to a difficult optimization problem that depends on the planar set, , supporting a peculiar positive function: any two
By assumption, the Radon transform depends only on , so its integral with respect to is equal to the value of the Radon transform at angle times the length of the support interval. Consequently, if the postulated function exists, then its Radon transform is otherwise (30) to integrate to one. where we have arbitrarily normalized This normalization can be done without loss of generality: if integrated to zero, then the Radon transform would be . identically zero, and so would We apply the inverse Radon transform formula [5] to (30) to obtain an expression for the postulated function in (31),
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
(a)
(b)
1047
(c)
Fig. 8. Assumed Radon transform for two angles, = 5=6, and = =2 implies that the function must be equal is equal to zero in the region shown in (c). The triangle cannot support a function having the special Radon transform property.
shown at the bottom of the page, where the asterisk denotes in the convolution with respect to . For all values of interior of , both bracketed quantities in the denominator of (31) are positive, and the function is strictly positive. When is outside of , we have to integrate through the singularities, and it is not obvious what happens. In particular, is a square, one can carry out a somewhat tedious if integration to show that the function does not vanish outside of (in fact, the function is negative in places), a proof by contradiction that the square does not support the postulated function.
(a)
(b)
(c)
(d)
B. Geometrical Considerations For elementary sets, simple geometric arguments demonstrate that the postulated function cannot exist. Recall again integrates to one. Suppose that the postulated function that exists. Then (30) implies that the integral over the region between any parallel lines is equal to the distance between the lines, divided by the length of the corresponding support interval. For the case of the triangle in Fig. 8, the assumed Radon , , gives the transform for the two angles, integral of the function within the subsets that are shown in Fig. 8(a) and (b). Together these facts imply that the function is equal to zero in the region shown in Fig. 8(c). A line integral through this zero-valued region is equal to zero, which contradicts the properties of the postulated function. In other words one can find three angles such that the assumed Radon transform (30) leads to a contradiction. For the case of a square, Fig. 9 shows that the assumed Radon transform is self-contradictory for four angles. C. Constructive Proof That
(e) Fig. 9. Assumed Radon transform for two angles: (a) = 3=4 and (b) = =4, implies (c), which when combined with (d), the assumed Radon transform for = 0, implies (e). The square cannot support a function having the special Radon transform property.
moment with respect to
of both of these expressions:
(32) For
, (32) implies, after some simplification, that
(33)
Is an Ellipse
be some convex, compact set, and assume that it Let supports a function having the Radon transform (30). We will must be elliptical. prove that We have the definition of the Radon transform (29) as well as the assumed Radon transform (30). We take the th
We denote the centroid of the function by that the function integrates to one):
(recalling
(34)
(31)
1048
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Let be the set obtained by shifting the set, , so that : the centroid is located at the origin, . Then the support function for is (35) and for Now solving (35) for (33), and using (34) gives
, substituting into (36)
So the set, , is symmetric about the origin, and the set, , must be symmetric about the centroid. This result alone implies that triangles and pentagons, for example, cannot support the postulated function. , (32) implies that For
3-D matched filters can reduce the size of a filterbank by a factor of five. For the elliptical velocity set we exploited a special property of the ellipse to obtain an analytical solution for the minimax filter: a convex elliptical region of the plane supports a positive-valued function such that any parallel lines that intersect the ellipse generate equal line integrals. While attempting to extend this solution technique to nonelliptical sets, we obtained the remarkable result that this property of the elliptical set is unique. No other planar convex, compact set supports a nonzero-valued function whose Radon transform depends only on the slope of the integrating line. ACKNOWLEDGMENT The authors would like to thank R. Burridge, P. C. Fishburn, L. Flatto, and T. J. Richardson for their helpful comments during this research.
(37) REFERENCES and in terms of Again expressing substituting into (37), using (34) and (36), gives
and
(38) is symmetric about the origin, its support function is nonnegative. Then (38) can be rewritten in terms of a quadratic form
Since
where the 2 2 matrix is symmetric, nonnegative definite. has the same After performing a coordinate rotation, form as the support function for the ellipse, (19). Consequently must be an elliptical set. V. CONCLUSIONS Detecting a moving object in a noisy image sequence requires a 3-D filtering operation to improve the SNR. If the object is not resolved by the sensor, then the effective object intensity distribution is the known point-spread function. However, in most applications the velocity of the object is unknown, and one needs a filter that works well over a set of velocities rather than for just one velocity. Under these conditions we have shown that an optimized filter, whose minimum SNR for a specified set of velocities is maximized, can be significantly more efficient than a 3-D matched filter: a SNR improvement of 3.5 dB for a circular disk, and for an elliptical region, an unbounded improvement as the eccentricity approaches one. For an elliptical velocity set (and its special cases, the disk and the line segment) we have a closed form analytical solution for the minimax filter. For other velocity sets we have obtained a useful constrained minimax filter, in closed form, whose SNR can be within 0.4 dB of the minimax SNR. Using constrained minimax filters rather than
[1] Y. Barniv, “Application of velocity filtering to optical-flow passive ranging,” IEEE Trans. Aerosp. Electron. Syst., vol. 28, pp. 957–969, 1989. [2] M. Born and E. Wolf, Principle of Optics, 6th ed. White Plains, NY: Pergamon, 1986. [3] L. T. Bruton and N. R. Bartley, “The enhancement and tracking of moving objects in digital images using adaptive three-dimensional recursive filters,” IEEE Trans. Circuits Syst., vol. 33, pp. 604–612, 1986. [4] Y. Chen, “On suboptimal detection of 3-dimensional moving targets,” IEEE Trans. Aerosp. Electron. Syst., vol. 25, pp. 343–350, 1989. [5] S. R. Deans, The Radon Transform and Some of Its Applications. New York: Wiley, 1983. [6] P. Embree, J. P. Burg, and M. M. Backus, “Wide-band velocity filtering—The pie-slice process,” Geophysics, vol. 28, pp. 948–974, 1963. [7] J. P. Fail and G. Grau, “Les filtres en eventail,” Geophysic. Prospect., vol. 11, pp. 131–163, 1963. [8] H. W. Guggenheimer, Differential Geometry. New York: McGrawHill, 1963. [9] K. Hsu and T. L. Marzetta, “Velocity filtering of acoustic well logging waveforms,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 37, pp. 265–274, 1989. [10] D. G. Luenberger, Linear and Nonlinear Programming, 2nd ed. Reading, MA: Addison-Wesley, 1984. [11] T. L. Marzetta, “Uniformly optimal 3-D fan filters for optical moving target detection,” in 1993 IEEE Int. Conf. Acoust., Speech, Signal Processing, Minneapolis, MN, Apr. 1993, pp. 543–545. [12] , “Fan filters, the 3-D Radon transform, and image sequence analysis,” IEEE Trans. Image Processing, vol. 3, pp. 253–264, May 1994. [13] , “Optimal detection of known moving objects in a noisy image sequence with velocity uncertainty,” in 1994 IEEE Int. Conf. Image Processing, Austin, TX, Nov. 1994, pp. 353–357. [14] S. C. Pohlig, “An algorithm for detection of moving optical targets,” IEEE Trans. Aerosp. Electron. Syst., vol. 25, pp. 56–63, 1989. [15] B. Porat and B. Friedlander, “A frequency domain algorithm for multiframe detection and estimation of dim targets,” IEEE Trans. Pattern Anal. Machine Intell., vol. 12, pp. 398–401, 1990. [16] I. S. Reed, R. M. Gagliardi, and H. M. Shao, “Applications of three dimensional filtering to moving target detection,” IEEE Trans. Aerosp. Electron. Syst., vol. 19, pp. 898–905, 1983. [17] I. S. Reed, R. M. Gagliardi, and L. Stotts, “Optical moving target detection with 3-D matched filtering,” IEEE Trans. Aerosp. Electron. Syst., vol. 24, pp. 327–335, 1988. [18] J. N. Sanders, “A method of determining filter spacing in assumed velocity filterbanks,” IEEE Trans. Aerosp. Electron. Syst., vol. 29, pp. 1043–1046, 1993. [19] G. Strang, Introduction to Applied Mathematics. Wellesley, MA: Wellesley-Cambridge Press, 1986. [20] H. L. VanTrees, Detection, Estimation, and Modulation Theory—Part I. New York: Wiley, 1968.
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
Thomas L. Marzetta (M’77–SM’93) received the Ph.D. degree in electrical engineering from the Massachusetts Institute of Technology, Cambridge. He is with the Mathematical Sciences Research Center, Bell Laboratories, Lucent Technologies, Murray Hill, NJ, where his principle focus is high capacity wireless communications. Prior to joining Bell Laboratories, he worked in petroleum exploration at Schlumberger-Doll Research, Ridgefield, CT, and in aerospace/electronics at Nichols Research Corporation, Wakefield, MA. Dr. Marzetta has been active in the IEEE Signal Processing Society as a Charter Member of the Technical Committee on Multidimensional Signal Processing, an Associate Editor of the TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, an Associate Editor of the TRANSACTIONS ON IMAGE PROCESSING, and as a Member of the Signal Processing Steering Committee of the IEEE Central New England Council. He was the recipient of the 1981 ASSP Paper Award.
1049
Larry A. Shepp was born in Brooklyn, NY. He received the Ph.D. degree from Princeton University, Princeton, NJ, in 1961, under W. Feller, with a dissertation on random walks. He is a Professor of Statistics at Rutgers University, New Brunswick, NJ. He worked at Bell Laboratories, Murray Hill, NJ, from 1962 to 1996 and spent sabbaticals in many places including Australia, Russia, Venezuela, MIT, and Berkeley. He is active in medical imaging, probability theory, mathematical finance, and combinatorics. Dr. Shepp is a member of NAS, IOM, IMS, and AMA.
1039
A Surprising Radon Transform Result and Its Application to Motion Detection Thomas L. Marzetta, Senior Member, IEEE, and Larry A. Shepp
Abstract—An elliptical region of the plane supports a positivevalued function whose Radon transform depends only on the slope of the integrating line. Any two parallel lines that intersect the ellipse generate equal line integrals of the function. We prove that this peculiar property is unique to the ellipse; no other convex, compact region of the plane supports a nonzero-valued function whose Radon transform depends only on slope. We motivate this problem by considering the detection of a constant-velocity moving object in a sequence of images, in the presence of additive, white, Gaussian noise. The intensity distribution of the object is known, but the velocity is only assumed to lie in some known set, for example, an ellipse or a rectangle. The object is to find a space-time linear filter, operating on the image sequence, whose minimum output signalto-noise ratio (SNR) for any velocity in the set is maximized. For an ellipse (and its special cases, the disk and the linesegment) the special Radon transform property of the ellipse enables us to obtain a closed-form, analytical solution for the minimax filter, which significantly outperforms the conventional three-dimensional (3-D) matched filter. This analytical solution also suggests a constrained minimax filter for other velocity sets, obtainable in closed form, whose SNR can be very close to the minimax SNR. Index Terms—Convex set, motion detection, Radon transform.
I. INTRODUCTION
A
CHALLENGING problem in image sequence analysis is to detect a moving object in the presence of additive noise [1], [3], [4], [14]–[18]. It is convenient and often reasonable to make several assumptions: that the object is moving with constant velocity, that its image intensity distribution does not change with time, and that the noise is zero mean Gaussian that is uncorrelated over both space and time. We denote the where represents time noisy image sequence by is Cartesian position within the image frame and (1) is the intensity distribution of the object, is the velocity of the object expressed in some consistent set of units, for example meters/s or pixels/frame, is the position of the object at reference time , and is additive white noise. By assumption, the noise variance is so large that the object cannot be detected
where
Manuscript received June 24, 1997; revised April 3, 1998. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Steven D. Blostein. T. L. Marzetta is with Bell Laboratories, Lucent Technologies, Murray Hill, NJ 07974-0636 USA (e-mail: [email protected]). L. A. Shepp is with the Department of Statistics, Rutgers University, New Brunswick, NJ 08903 USA (e-mail: [email protected]). Publisher Item Identifier S 1057-7149(99)06116-3.
reliably merely by thresholding the image sequence, or by applying a two-dimensional (2-D) (spatial) matched filter to each image frame. In general we need a more sophisticated three-dimensional (3-D) detection scheme that operates over multiple pixels and multiple frames of the image sequence. A. Four Motion Detection Problems The nature of the detection problem depends on the extent of prior information that is available concerning the object of interest. In almost any practical problem the reference position , is unknown. The velocity may or may of the object, not be known, and the object intensity distribution may or may not be known. It is instructive to consider briefly the four cases, each of which requires a distinct detection strategy. The fourth case, an unknown velocity and a known object intensity distribution, is the most interesting case, as well as the focus of this paper. 1) Velocity Known; Object Intensity Distribution Known: Here we know everything about the signal except for the reference position. Correlating the noisy image sequence with a replica of the signal to be detected for all possible generates a set of sufficient statistics. The actual detection operation compares the largest value of the correlation with a threshold; this constitutes a generalized likelihood ratio test [20] for the presence or absence of the moving object. It is necessary to limit the temporal duration of the correlation operation if only because the object remains in the field of view for a finite length of time. The 3-D (i.e., space-time) correlation operation is equivalent to applying the so-called 3-D matched filter [16], whose impulse response is (2) is a finite duration temporal window, for example where a rectangle or a raised cosine. When applied to the noisy image sequence the filter produces an image sequence containing a moving object having the original velocity, and with an improved signal-to-noise ratio (SNR). The 3-D matched filter operates in two ways: it implements a spatial matched filter, and it integrates the image sequence along lines in the space that coincide with the trajectory of the moving object. The signal-to-noise improvement is proportional to the duration of the temporal window. 2) Velocity Known; Object Intensity Distribution Unknown: is unknown, and that the image sequence Suppose that is observed over the time interval, [0, ]. It is straightforward to show that, with knowledge only of the velocity, a sufficient
1057–7149/99$10.00 1999 IEEE
1040
Fig. 1. Bank of 3-D matched filters can cover a region of possible object velocities. Filters achieve their maximum SNR at the grid points; velocity mismatch occurs between grid points.
statistic is generated by integrating the image sequence along the trajectory of the moving object. This sufficient statistic is a single image that contains the object with an enhanced SNR that is proportional to
Therefore the original detection problem reduces to detecting an unknown object in a single noisy image. To proceed any . further it is necessary to assume something about For example if the extent of the object were known, one with a disk or rectangle whose size could correlate , is comparable to that of the object. The generation of and the spatial correlation operation together are equivalent to is replaced applying the 3-D matched filter (2), where by a disk or rectangle within the formula. 3) Velocity Unknown; Object Intensity Distribution Unknown: Situations where the velocity is unknown occur more frequently than where the velocity is known. First consider and the velocity are unknown. the case where both Conceptually one could repeatedly apply the 3-D matched filter for all possible velocities—an impractical scheme. In practice a bank of 3-D matched filters has been proposed as in Fig. 1, whose design velocities span the possible range of velocity of the object [1], [15], [17]. Given the severe computational burden of implementing a 3-D filter it is desirable to maximize the efficiency of the filterbank to achieve reliable detection with as few filters as possible by optimizing the temporal window and the design velocities. When the object velocity falls between velocity-grid points in the filterbank velocity mismatch occurs, and the 3-D matched filter is no longer optimal; in fact, increasing the temporal duration of the filter indefinitely eventually results in a reduced SNR, because then the filter integrates mostly noise. The optimization of a bank of 3-D matched filters can be difficult because there is no general closed-form expression for SNR improvement when there is velocity mismatch. An alternative approach handles velocity uncertainty by designing the filter to have an adequate SNR over a set of velocities—for example, a square or a hexagon—rather than
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Fig. 2. Filterbank where each filter yields the stipulated SNR uniformly for every velocity inside of its tile.
to have a high SNR at merely one velocity [11]–[13]. One can design a bank of such filters for a larger velocity set by tiling the region of interest in the velocity plane as in Fig. 2, such that each tile corresponds to one of the filters. For the case where the velocity is unknown and the object intensity distribution is unknown, one can require that a 3-D filter have minimum white-noise sensitivity subject to the constraint that it pass all moving objects, whose velocity is (that corresponds to one of the tiles contained in a set in Fig. 2) with unity gain [11], [12]. This design procedure yields a constant SNR over the velocity set. The resulting filter is an example of a fan filter—originally used in geophysical applications to filter wavefields on the basis of plane wave slowness [6], [7], [9]. Since the fan filter has unity gain over the 3-D frequency support of the moving object, it does not destroy any information. For elementary velocity sets there are closed-form expressions for the frequency response of the filter and for its SNR improvement. 4) Velocity Unknown; Object Intensity Distribution Known: In certain cases, the object intensity distribution is known, but the velocity of the object is unknown. This can occur, for example, when the angular extent of the object is much smaller than the angular resolution of the sensor, in which case the effective object intensity distribution is equal to the known point spread function of the sensor. An intuitively appealing approach to this problem would cascade a spatial matched filter (e.g., matched to the known object intensity distribution) with the 3-D fan filter mentioned in Section I-A3. It turns out, however, that this combination can be improved upon. Instead, without imposing any a priori structure one can formulate a minimax design problem for the filter where, for all velocities in the set , the minimum SNR is maximized [13]. A general analytical solution to the minimax problem was not found, but for the case where is an elliptical region a closed-form analytical solution was found whose existence depends on the set supporting a peculiar function. B. Peculiar Radon Transform Property In the case of an elliptical velocity set, the key to the analytical solution for the minimax filter is the existence of a positive Lagrange multiplier function on the elliptical region whose Radon transform depends only on the slope of the
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
1041
sequence yields the best SNR in a minimax sense for all velocities in the set, . The noise variance at the output of the . The result of filter is equal to filtering the moving object is again a moving object that has the same velocity, but a different object intensity distribution. We want to maximize the amplitude of the filtered moving object, which is proportional to the filter gain, at the arbitrary point, . (Shifting this point would merely shift the optimized filter impulse response.) The SNR at the output of the filter is a function of the unknown velocity, SNR Fig. 3. Convex elliptical region supports a positive-valued function such that all parallel lines that intersect the ellipse generate equal line integrals; the Radon transform depends only on the slope of the integrating line.
integrating line; any two parallel lines that intersect the ellipse generate equal line integrals of the function (see Fig. 3). It is somewhat surprising that such a function exists. It is even more surprising, as we show in this paper, that the elliptical region is the only convex compact set that supports such a function. C. Organization of Paper Section II formulates the design of the minimax filter, which is equivalent to an infinite dimensional quadratic programming problem. We apply the Kuhn–Tucker conditions to obtain the optimal filter in terms of a nonnegative-valued Lagrange , defined on the velocity set multiplier function, which, in turn, satisfies a 2-D linear integral equation. While a general solution to the integral equation eludes us, we show that an analytical solution can be obtained if a function exists that has the peculiar Radon transform property mentioned in Section I-B. There is such a function when is an ellipse. The minimax filter for the ellipse suggests a constrained minimax filter for other sets, such as the rectangle or the hexagon that, although suboptimal, can be obtained analytically, and whose performance approaches that of the optimal minimax filter. Section III compares the performance of the minimax filter, the constrained minimax filter, and the optimized 3-D matched filter for some cases of interest. Section IV provides a constructive proof that the ellipse is the only convex, compact set that supports a nonzero-valued function whose Radon transform depends only on the slope of the integrating line. II. MINIMAX 3-D FILTER
(3) where the numerator is the square of the filter gain as a function of velocity. The optimal filter is the solution to the following minimax optimization problem: SNR This minimax problem is equivalent to an infinite dimensional quadratic programming problem where we minimize the noise variance at the output of the filter subject to a linear inequality constraint on the filter gain, for all possible velocities, subject to
(4) where is an arbitrary constant (changing merely scales the amplitude of the filter without changing the SNR). It is convenient to reformulate the minimax problem in the 3-D frequency domain. The 3-D frequency response of the filter is
where we use the geophysical convention with respect to the signs in the exponent. The 3-D Fourier transform of the moving object is
A. Problem Statement We observe the noisy image sequence, (1), where the object , is known, the reference position intensity distribution, , is unknown, and the velocity of the of the object, , is unknown. However the velocity is assumed object, to lie in some known convex compact set, . The noise, , is zero mean, Gaussian, white with known spectral . density, We seek a 3-D linear shift-invariant filter having impulse , which when applied to the noisy image response,
where is the Dirac delta function, and is the 2-D Fourier transform of the object intensity distribution
1042
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Thus, the moving object occupies a plane in the 3-D frequency space that passes through the origin
(7)
(5) With the above notation the minimax problem (4) becomes (8) subject to
(9)
(6) A feasible solution (e.g., a finite energy filter that satisfies the constraints) always exists: an example is the cascade of the spatial matched filter and the fan filter mentioned in Section I-A4. For mathematical convenience, we have formulated the detection problem under the assumption that the noisy image “sequence” is measured continuously with respect to space and time, and the 3-D filter is a continuous parameter filter rather than a sampled data filter. In the case of a diffraction-limited optical sensor, the modulation transfer function is spatially bandlimited, and so is the object intensity distribution [2]; it follows from the relation (5) that for all finite velocities, the moving object is temporally bandlimited as well. Under these conditions a spatially and temporally lowpass filtered version of the continuous parameter noisy image sequence could be sampled at the Nyquist rate without loss of information, and the discrete filter would merely be a sampled version of the continuous filter. B. Kuhn–Tucker Conditions The inequality constraints involved in the optimization (6) can be converted into equality constraints by subtracting a nonnegative function of velocity from the left-hand side of the constraint; in turn this nonnegative function is equal to the . We include square of an unconstrained function, the resulting equality constraints in the minimization through the method of Lagrange multipliers. Since there are an infinite number of constraints we use a Lagrange multiplier function, . The resulting minimization problem with respect , , and becomes to
The calculus of variations gives three conditions:
The first condition (7) gives the optimal filter in terms of the Lagrange multiplier function, and the second condition (8) is merely the equality constraint for every velocity. The substitution of (7) into (8), and the integration over , yields a 2-D integral equation for the Lagrange multiplier function
(10) appears in the integral equation, Note that the function and are also mutually dependent through the and that third condition (9). The quantitative interpretation of the Lagrange multiplier is that it is the derivative of the optimized cost (e.g., the noise variance at the output of the filter) with respect to the righthand side of the constraint, [19]. As increases, implying an increase in filter gain, the noise variance at the output of the optimized filter cannot decrease. Consequently, the Lagrange multiplier function must be nonnegative. We combine this fact with (9) to obtain the following Kuhn–Tucker conditions [10]:
if
then
if
then
(11) .
is equal to zero, the filter gain is For velocities where exactly equal to and the filter achieves the minimax SNR is not equal to zero, exactly, while for velocities where the filter gain exceeds and the SNR exceeds the minimax SNR. The implications of the Kuhn–Tucker conditions are the following. must be 1) The solution of the integral equation for nonnegative. is positive the filter yields 2) For all velocities where exactly its minimax SNR. 3) For all velocities where the SNR exceeds the minimax must be equal to zero. SNR, In short, constraints are active where the Lagrange multiplier function is greater than zero. We do not have an analytical solution to the integral equation in combination with the Kuhn–Tucker conditions for the general case. However, in the next two subsections we derive some general properties of the optimal filter, and we actually obtain a closed-form analytical solution for the case is an ellipse. where
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
1043
Then one can verify by direct substitution that the solutions in the case of the translated set are
Fig. 4. Support function for a convex, compact set; distance between the tangent line and the origin.
f ( )
is the signed
C. General Properties of Minimax Filter The minimax filter (7) is equal to the product of the spatial , and a 3-D filter that consists of a linear matched filter, superposition of impulsive planes that are defined by (5). As a result, the spatial frequency support of the minimax filter is identical to that of the object intensity distribution, and , the union of the impulsive planes for a particular determines the temporal frequency support. We denote the , temporal frequency support by for some In turn this set can be conveniently expressed in terms of the support function [8] for the convex velocity set, (12) As shown in Fig. 4, the support function is equal to the signed distance between a tangent line and the origin. It follows easily that (13) Thus the convex set is defined by a pair of parallel tangent lines for every angle . We express spatial frequency in polar coordinates, (14) It follows that the temporal frequency support for the minimax filter is the continuous interval
The shift of the Lagrange multiplier function changes the frequency response (7) of the filter nontrivially, however a direct calculation of the energy of the filter shows that its noise sensitivity, and therefore the minimax SNR, is unchanged. If one is designing a bank of minimax filters using identical tiles, then only one optimization problem needs to be solved because of this property, and all of the filters yield the same minimax SNR. Property 2—For a Fixed Shape and Orientation of the Velocity Set, the Minimax SNR Is Inversely Proportional to the and Square-Root of the Area of the Set: Suppose that satisfy the integral equation and the Kuhn–Tucker conditions for some velocity set, . Let be the velocity set obtained by expanding (in both dimensions) the set by the factor , which increases the area by the factor ,
It can be shown that the solutions in the case of the expanded are set
A calculation of the energy of the resulting minimax filter discloses that the minimax SNR (in units of power/power) is decreased by the factor . For example, to increase the minimax SNR by 3 dB (i.e., doubling the SNR), we have to decrease the linear size of the velocity set by a factor of two. As a result the velocity-area of each tile is reduced by a factor of four, or equivalently the number of filters has to increase by a factor of four to cover the same region of interest. D. Closed-Form Solution When
Is Elliptical
In general we can replace the integral equation involving by an equivalent integral equation involving the Radon . The inner integral of (10) is related to the transform of Radon transform of the Lagrange multiplier function; with the polar representation of spatial frequency, (14) we have
(15) Without actually solving for the minimax filter, we can deduce two important properties that greatly simplify the task of designing a filterbank. Property 1—The Minimax SNR Is Invariant to Translations of the Velocity Set: Suppose that for the set , we have and that satisfy the integral equation (10) solutions for and the Kuhn–Tucker conditions (11). Let be the velocity set , i.e., obtained by translating the set by the amount
where
is the Radon transform of
1044
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
After substituting this expression into (10) the integral equation becomes
The calculated SNR is consistent with Property 2 of Section IIC. We show in Section IV that the ellipse is the only set that has the peculiar Radon transform property that we exploited here to obtain an analytical solution for the minimax filter.
(16) Note that the velocity affects the integral only through the argument of the Radon transform. If we can find a nonnegative Lagrange multiplier function whose Radon transform has no -dependence, for all values of such that the integrating line intersects the velocity set, then a scaled version of automatically satisfies (16) with . For example suppose that is the elliptical region, . Then it can be verified by elementary calculation that the following positive function:
E. Constrained Minimax Filter for Arbitrary Convex Velocity Sets
(17) has the Radon transform
For cases where the velocity set is not an ellipse, one could always obtain the minimax filter by numerical means. To that end, a dual optimization problem is more tractable than the primal optimization problem (6). It can be shown that the dual problem involves choosing the nonnegative-valued Lagrange multiplier function to maximize a convex quadratic functional, with no other constraints. However we do not pursue the numerical approach in this paper. Instead, we obtain a useful suboptimal filter analytically. For the case where is an ellipse, note that, within its passband, the frequency response of the minimax filter (20) is independent of . This motivates us to consider, in general, a filter of the form
otherwise
otherwise (18) where
is the support function, (12), (19)
depends only on spatial frequency, and where where again spatial frequency is expressed in polar coordinates. The SNR of this filter is the same for all velocities within . The that maximizes the SNR gives (to within a choice of multiplicative constant)
After substituting (17) into (7) and using (18), we obtain the closed form expression for the minimax filter (to within an arbitrary multiplicative constant)
(22) otherwise. We refer to this as the constrained minimax filter. The resulting SNR is
otherwise. (20)
SNR
The frequency support of the filter is an elliptical cone. The is no cause for alarm. In the singularity at , vicinity of the origin the factor, , is both absolutely- and square-integrable with . The presence of this singularity is inturespect to itively reasonable: the frequency support plane (5) associated with every velocity in passes through the origin, so it makes sense to concentrate the filter energy near the origin rather than to spread the filter energy evenly throughout the filter support. The resulting SNR is the same for all velocities inside of the ellipse:
(23) It is easily shown that the two properties of Section II-C apply to the constrained minimax filter as well. We later show that the performance of this filter can be very close to that of the optimal minimax filter. III. PERFORMANCE RESULTS In what follows we assume that the object intensity distribution is Gaussian-shaped:
SNR (21)
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
1045
(24) is the integrated intensity of the where object, and represents the size of the object. This is a common model for the point-spread function of an optical sensor. Additionally it turns out that this object intensity distribution allows us to obtain performance results analytically. A. Disk-Shaped Velocity Set: Minimax Filter Versus 3-D Matched Filter Here we assume a disk-shaped velocity set of radius centered at the origin, and we compare the performance of the minimax filter with that of an optimized 3-D matched filter. 1) Minimax Filter: The minimax filter is given by (20), , and the minimax SNR is given evaluated for by (21), easily evaluated analytically for the Gaussian object intensity distribution (24). We can also evaluate the SNR for velocities outside of the disk, which is less than the minimax SNR, SNR
. (25) 2) 3-D Matched Filter: The 3-D matched filter (2) is characterized by its design velocity and by its temporal window. Its 3-D frequency response, when tuned for zero-velocity (the center of ), is
Because of the circular symmetry of , the resulting SNR depends only on the magnitude of the velocity, and it decreases monotonically with the speed. We optimize the filter in a to maximize the SNR on minimax sense by choosing the boundary of the velocity disk [specifically at the point, ]. The result of the optimization is
Fig. 5. Filterbank based on hexagonal tiling. Each tile is covered by a constrained minimax filter that provides uniformly good SNR. The SNR of the filters is within 0.4 dB of the minimax SNR.
improvement factor is (0.5 dB).1 Equivalently, for the same guaranteed SNR, the minimax filter can cover a velocity disk whose radius is 2.2 times greater than that of the 3-D matched filter. The improved performance of the minimax filter relative to the 3-D matched filter is even more dramatic for an ellipse. The 3-D matched filter has to be designed for the major axis , and its minimum SNR within the of the ellipse, say ellipse is inversely proportional to . In contrast, the minimax remains SNR increases without limit as decreases while constant. B. Designing a Filterbank: Constrained Minimax Filters versus 3-D Matched Filters Suppose that the possible range of velocity of the object is so great that the required SNR cannot be realized by a single filter, necessitating the design of a bank of filters. We compare two approaches. The first approach tiles the velocity plane with hexagons, with each tile covered by the constrained minimax filter (22). The second approach uses a bank of optimized 3-D matched filters whose design velocities are on a hexagonal grid. 1) Constrained Minimax Filters for Hexagonal Tiles: Fig. 5 shows a hexagonal tiling of the velocity plane, with the halfspacing between the centers of adjacent tiles equal to . The support function for a hexagon centered at the origin is mod . An evaluation of (23) with this support function gives the resulting SNR for any of the constrained minimax filters when the velocity is in the corresponding tile; this expression is then solved for the required tile size:
The SNR resulting from the optimized 3-D matched filter is SNR
SNR or
(26) 3) Comparison: Comparing the SNR’s of the two filters, (25) and (26), we find that the minimax filter performs better than the 3-D matched filter everywhere in the velocity disk: on the boundary of the disk the SNR improvement factor is (3.5 dB), while at the center of the disk the
SNR
(27)
It would be interesting to compare the SNR for the constrained minimax filter with the SNR for the ideal minimax 1 An even greater improvement factor was given in [13], however a different object intensity distribution than (24) was assumed, and the 3-D matched filter was not truly optimized.
1046
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Fig. 6. Performance of minimax filter for hexagon is bounded by that of minimax filters for outer disk and inner disk.
filter, but this would require a numerical optimization. Instead we can obtain upper and lower bounds on the minimax SNR, equal to the minimax SNR’s for the two disks shown in Fig. 6, , and the inner disk the outer disk having a radius of having a radius of : SNR A comparison of the upper bound with (27) shows that the optimal minimax filter for the hexagon would have an SNR (0.42 dB) that is at most a factor of greater than that of the constrained minimax filter. The ideal minimax filter could reduce the size of the filterbank by, at most, 20%. 2) 3-D Matched Filters on a Hexagonal Velocity Grid: Fig. 7 illustrates a bank of 3-D matched filters arranged on a hexagonal velocity grid. The half-spacing between adjacent grid points is . Each filter has a disk-shaped region of , in which the filter achieves coverage, of radius the stipulated SNR. Of necessity, the disks have to be overlapped in the manner shown. Equation (26) relates to the minimum SNR that is achieved by the filterbank; this expression, in turn, gives the minimum SNR in terms of , SNR
Fig. 7. Bank of 3-D matched filters is less efficient than bank of constrained minimax filters shown in Fig. 5. Each filter achieves the stipulated SNR on the boundary of the disk, with maximum SNR at the center of the disk. Disks must overlap to compensate for velocity mismatch.
parallel lines that intersect the set generate equal line integrals of the function. We now prove that this property is unique to the elliptical region. A. Desired Radon Transform Property be a convex, compact set in the plane, and let be a function defined on . The Radon transform of the function is the line integral of the function, for the , integrating line, Let
(29) The Radon transform vanishes for outside of the support in, where is the support function, terval, (12). , whose The object is to find a nonzero-valued function, for all values of Radon transform is independent of within the support interval. If such a function exists, then it is straightforward to infer its Radon transform. To do so, we integrate both sides of (29) with respect to , as follows:
or SNR
(28)
3) Comparison: Upon comparing (27) with (28) we find that the constrained minimax filterbank is considerably more efficient than the 3-D matched filterbank. To achieve the same minimum SNR with 3-D matched filters we must have
Using the 3-D matched filter rather than the constrained times as many minimax filter would require filters. IV. NEW RADON TRANSFORM RESULT In the previous section, we found a closed-form analytical solution to a difficult optimization problem that depends on the planar set, , supporting a peculiar positive function: any two
By assumption, the Radon transform depends only on , so its integral with respect to is equal to the value of the Radon transform at angle times the length of the support interval. Consequently, if the postulated function exists, then its Radon transform is otherwise (30) to integrate to one. where we have arbitrarily normalized This normalization can be done without loss of generality: if integrated to zero, then the Radon transform would be . identically zero, and so would We apply the inverse Radon transform formula [5] to (30) to obtain an expression for the postulated function in (31),
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
(a)
(b)
1047
(c)
Fig. 8. Assumed Radon transform for two angles, = 5=6, and = =2 implies that the function must be equal is equal to zero in the region shown in (c). The triangle cannot support a function having the special Radon transform property.
shown at the bottom of the page, where the asterisk denotes in the convolution with respect to . For all values of interior of , both bracketed quantities in the denominator of (31) are positive, and the function is strictly positive. When is outside of , we have to integrate through the singularities, and it is not obvious what happens. In particular, is a square, one can carry out a somewhat tedious if integration to show that the function does not vanish outside of (in fact, the function is negative in places), a proof by contradiction that the square does not support the postulated function.
(a)
(b)
(c)
(d)
B. Geometrical Considerations For elementary sets, simple geometric arguments demonstrate that the postulated function cannot exist. Recall again integrates to one. Suppose that the postulated function that exists. Then (30) implies that the integral over the region between any parallel lines is equal to the distance between the lines, divided by the length of the corresponding support interval. For the case of the triangle in Fig. 8, the assumed Radon , , gives the transform for the two angles, integral of the function within the subsets that are shown in Fig. 8(a) and (b). Together these facts imply that the function is equal to zero in the region shown in Fig. 8(c). A line integral through this zero-valued region is equal to zero, which contradicts the properties of the postulated function. In other words one can find three angles such that the assumed Radon transform (30) leads to a contradiction. For the case of a square, Fig. 9 shows that the assumed Radon transform is self-contradictory for four angles. C. Constructive Proof That
(e) Fig. 9. Assumed Radon transform for two angles: (a) = 3=4 and (b) = =4, implies (c), which when combined with (d), the assumed Radon transform for = 0, implies (e). The square cannot support a function having the special Radon transform property.
moment with respect to
of both of these expressions:
(32) For
, (32) implies, after some simplification, that
(33)
Is an Ellipse
be some convex, compact set, and assume that it Let supports a function having the Radon transform (30). We will must be elliptical. prove that We have the definition of the Radon transform (29) as well as the assumed Radon transform (30). We take the th
We denote the centroid of the function by that the function integrates to one):
(recalling
(34)
(31)
1048
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 8, AUGUST 1999
Let be the set obtained by shifting the set, , so that : the centroid is located at the origin, . Then the support function for is (35) and for Now solving (35) for (33), and using (34) gives
, substituting into (36)
So the set, , is symmetric about the origin, and the set, , must be symmetric about the centroid. This result alone implies that triangles and pentagons, for example, cannot support the postulated function. , (32) implies that For
3-D matched filters can reduce the size of a filterbank by a factor of five. For the elliptical velocity set we exploited a special property of the ellipse to obtain an analytical solution for the minimax filter: a convex elliptical region of the plane supports a positive-valued function such that any parallel lines that intersect the ellipse generate equal line integrals. While attempting to extend this solution technique to nonelliptical sets, we obtained the remarkable result that this property of the elliptical set is unique. No other planar convex, compact set supports a nonzero-valued function whose Radon transform depends only on the slope of the integrating line. ACKNOWLEDGMENT The authors would like to thank R. Burridge, P. C. Fishburn, L. Flatto, and T. J. Richardson for their helpful comments during this research.
(37) REFERENCES and in terms of Again expressing substituting into (37), using (34) and (36), gives
and
(38) is symmetric about the origin, its support function is nonnegative. Then (38) can be rewritten in terms of a quadratic form
Since
where the 2 2 matrix is symmetric, nonnegative definite. has the same After performing a coordinate rotation, form as the support function for the ellipse, (19). Consequently must be an elliptical set. V. CONCLUSIONS Detecting a moving object in a noisy image sequence requires a 3-D filtering operation to improve the SNR. If the object is not resolved by the sensor, then the effective object intensity distribution is the known point-spread function. However, in most applications the velocity of the object is unknown, and one needs a filter that works well over a set of velocities rather than for just one velocity. Under these conditions we have shown that an optimized filter, whose minimum SNR for a specified set of velocities is maximized, can be significantly more efficient than a 3-D matched filter: a SNR improvement of 3.5 dB for a circular disk, and for an elliptical region, an unbounded improvement as the eccentricity approaches one. For an elliptical velocity set (and its special cases, the disk and the line segment) we have a closed form analytical solution for the minimax filter. For other velocity sets we have obtained a useful constrained minimax filter, in closed form, whose SNR can be within 0.4 dB of the minimax SNR. Using constrained minimax filters rather than
[1] Y. Barniv, “Application of velocity filtering to optical-flow passive ranging,” IEEE Trans. Aerosp. Electron. Syst., vol. 28, pp. 957–969, 1989. [2] M. Born and E. Wolf, Principle of Optics, 6th ed. White Plains, NY: Pergamon, 1986. [3] L. T. Bruton and N. R. Bartley, “The enhancement and tracking of moving objects in digital images using adaptive three-dimensional recursive filters,” IEEE Trans. Circuits Syst., vol. 33, pp. 604–612, 1986. [4] Y. Chen, “On suboptimal detection of 3-dimensional moving targets,” IEEE Trans. Aerosp. Electron. Syst., vol. 25, pp. 343–350, 1989. [5] S. R. Deans, The Radon Transform and Some of Its Applications. New York: Wiley, 1983. [6] P. Embree, J. P. Burg, and M. M. Backus, “Wide-band velocity filtering—The pie-slice process,” Geophysics, vol. 28, pp. 948–974, 1963. [7] J. P. Fail and G. Grau, “Les filtres en eventail,” Geophysic. Prospect., vol. 11, pp. 131–163, 1963. [8] H. W. Guggenheimer, Differential Geometry. New York: McGrawHill, 1963. [9] K. Hsu and T. L. Marzetta, “Velocity filtering of acoustic well logging waveforms,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 37, pp. 265–274, 1989. [10] D. G. Luenberger, Linear and Nonlinear Programming, 2nd ed. Reading, MA: Addison-Wesley, 1984. [11] T. L. Marzetta, “Uniformly optimal 3-D fan filters for optical moving target detection,” in 1993 IEEE Int. Conf. Acoust., Speech, Signal Processing, Minneapolis, MN, Apr. 1993, pp. 543–545. [12] , “Fan filters, the 3-D Radon transform, and image sequence analysis,” IEEE Trans. Image Processing, vol. 3, pp. 253–264, May 1994. [13] , “Optimal detection of known moving objects in a noisy image sequence with velocity uncertainty,” in 1994 IEEE Int. Conf. Image Processing, Austin, TX, Nov. 1994, pp. 353–357. [14] S. C. Pohlig, “An algorithm for detection of moving optical targets,” IEEE Trans. Aerosp. Electron. Syst., vol. 25, pp. 56–63, 1989. [15] B. Porat and B. Friedlander, “A frequency domain algorithm for multiframe detection and estimation of dim targets,” IEEE Trans. Pattern Anal. Machine Intell., vol. 12, pp. 398–401, 1990. [16] I. S. Reed, R. M. Gagliardi, and H. M. Shao, “Applications of three dimensional filtering to moving target detection,” IEEE Trans. Aerosp. Electron. Syst., vol. 19, pp. 898–905, 1983. [17] I. S. Reed, R. M. Gagliardi, and L. Stotts, “Optical moving target detection with 3-D matched filtering,” IEEE Trans. Aerosp. Electron. Syst., vol. 24, pp. 327–335, 1988. [18] J. N. Sanders, “A method of determining filter spacing in assumed velocity filterbanks,” IEEE Trans. Aerosp. Electron. Syst., vol. 29, pp. 1043–1046, 1993. [19] G. Strang, Introduction to Applied Mathematics. Wellesley, MA: Wellesley-Cambridge Press, 1986. [20] H. L. VanTrees, Detection, Estimation, and Modulation Theory—Part I. New York: Wiley, 1968.
MARZETTA AND SHEPP: SURPRISING RADON TRANSFORM RESULT
Thomas L. Marzetta (M’77–SM’93) received the Ph.D. degree in electrical engineering from the Massachusetts Institute of Technology, Cambridge. He is with the Mathematical Sciences Research Center, Bell Laboratories, Lucent Technologies, Murray Hill, NJ, where his principle focus is high capacity wireless communications. Prior to joining Bell Laboratories, he worked in petroleum exploration at Schlumberger-Doll Research, Ridgefield, CT, and in aerospace/electronics at Nichols Research Corporation, Wakefield, MA. Dr. Marzetta has been active in the IEEE Signal Processing Society as a Charter Member of the Technical Committee on Multidimensional Signal Processing, an Associate Editor of the TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, an Associate Editor of the TRANSACTIONS ON IMAGE PROCESSING, and as a Member of the Signal Processing Steering Committee of the IEEE Central New England Council. He was the recipient of the 1981 ASSP Paper Award.
1049
Larry A. Shepp was born in Brooklyn, NY. He received the Ph.D. degree from Princeton University, Princeton, NJ, in 1961, under W. Feller, with a dissertation on random walks. He is a Professor of Statistics at Rutgers University, New Brunswick, NJ. He worked at Bell Laboratories, Murray Hill, NJ, from 1962 to 1996 and spent sabbaticals in many places including Australia, Russia, Venezuela, MIT, and Berkeley. He is active in medical imaging, probability theory, mathematical finance, and combinatorics. Dr. Shepp is a member of NAS, IOM, IMS, and AMA.
