Binary Restoration of Thin Objects in Multidimensional ... - CiteSeerX
Binary Restoration of Thin Objects in Multidimensional Imagery Jeffrey E. Boyd Department of Electrical and Computer Engineering University of California, San Diego La Jolla, CA 92093-0407 [email protected] Jean Meloche Department of Statistics, University of British Columbia Vancouver, BC, Canada V6T 1Z2 [email protected]
Abstract
gado et al. [16]. Therefore, detection of concealed explosives poses a particularly difficult problem, reWe present a method for restoration of noisy tomo- quiring the ability to identify thin objects in images graphic images for detecting thin objects, such as ex- with high noise and limited resolution. plosives. Use of a weighted mean-square estimate The goal of the work presented here is to find a optimizes the solution to place emphasis on the infresolution to this problem. Thin objects in two dimenquent, but significant local structure associated with sions are lines or curves no more than one pixel wide, thin objects. Experimental results show successful and in three dimensions they are sheets no more than restoration at very high noise levels. one voxel thick. Ultimately, we wish to detect these thin objects in data where the standard deviation of the noise is one to two times the density of the ob1 Introduction jects of interest. An important feature of the underlying image is its binary nature. Tomographic images Detection of certain objects in multidimensional improvide spatial density maps, so that, within a linear age data is a common problem in industrial inspectransformation, the density of an object is one on a tion. The objects may be a normal part of the subject background of zero. The assumption that there are of the inspection, or they may be anomalies. Deteconly two possible density values that can be mapped tion is more difficult if the object has been deliberonto ones and zeros is significant and, if exploited ately concealed, which is likely with explosives or as shown here, leads to successful detection in very illegal drugs. A common method used to conceal noisy data. While this work focuses on detection of malleable material is to form it into a thin sheet. It is concealed explosives, the method has potential for well known that this is effective for conventional Xapplication in other inspection tasks, such as stress ray imagery where a thin object is almost invisible fractures, and disbonds and delaminations in comunless viewed on end. The use of tomography ofposite materials. fers a solution to the problem of imaging thin sheets by viewing the object from many angles, but tomoWe view the detection problem as one of image graphic systems sensitive to the presence of explo- restoration. That is, we assume that if the noisy data sives are subject to high noise levels and limited res- are restored to a nearly noise free state, then the deolution. For example, see Seed et al. [18] and Mor- tection problem becomes trivial. Such is the case
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
(a)
(b)
Figure 1: Comparison of wide and thin object detection: (a) signal (solid) and smoothed signal (dashed) for a wide object, and (b) signal and smoothed signal for thin object.
when the imaging modality is sensitive to some characteristic unique to the objects of interest [16, 18]. The simplest method of restoration is to smooth the data with an isotropic linear shift-invariant lowpass filter. Such an approach tends to blur discontinuities, which is unacceptable for detection of thin objects. Figure 1 illustrates this fundamental problem. Smoothing a wide object blurs its edges but the object is still clearly present. The thin object all but disappears after smoothing. Therefore, isotropic smoothing of the data is not viable. Another option is to use many non-isotropic low-pass filters, each tuned to smooth only in one orientation. The required number of filters depends on the spatial support of the filters. This method requires much more computation and leaves the problem of reconciling the output of the different filters. Perona [17] describes steerable filters that are linear combinations of several non-isotropic filters. The number of filters required for accurate results increases with the aspect ratio of the individual isotropic filters. In the case of thin objects, the aspect ratio is high and a large number of filters is required. Several methods of restoration consider the presence of discontinuities in order to avoid smoothing across them. Geman and Geman [5], in their seminal work, model an image as a Markov random field (MRF) and include edge processes to mark discontinuities and avoid blurring. Blake and Zisserman [1] proposed the graduated non-convexity (GNC) algorithm that approximates an image with discontinuities, later implemented in hardware by Harris et al. [10]. Geiger and Girosi [4] derive GNC from MRF’s. Leclerc [13] optimizes a restoration by minimizing
2
the length of a code that describes the image regions. Keeler [12] proposes a similar method that incorporates topological considerations. These methods all rely, at some point, on a difference technique to infer the presence of a discontinuity or edge. As a result, these methods are effective where noise levels are not too great, but when the noise levels are high, reliable detection of edges is nearly impossible and these methods fail. Mcdonald and Owen [14] and Hall and Titterington [8] avoid reliance on edge detection by using split filters. The method presumes the existence of a specific number of discontinuities in each neighborhood of the data, and then computes the optimum positions of discontinuities to minimize the error between the data and the between-discontinuity fits. In multiple dimensions, the filter optimizes by performing an exhaustive search over a large number of possible discontinuity orientations and positions, and is very time-consuming. We present a method for binary image restoration, based on Meloche and Zamar [15], to exploit the binary nature of some data and achieve superior restoration results. For any neighborhood in a binary image, there is only a finite number of possible configurations of ones and zeros. Binary image restoration encapsulates knowledge of the structure of the true image in the function Q, the frequency distribution over the set of possible neighborhood configurations. In effect, selection of Q tunes the restoration to the characteristics of the objects to be found. Binary restoration is applicable to the detection of thin sheets so long as the underlying image is binary in nature. Our results indicate that the method is effective in restoring thin objects in noisy data, and that it is not necessary to know the true Q, but it is better to use a weighted Q, based on prior information.
2 Binary Image Restoration The following is a summary of the binary image restoration method proposed by Meloche and Zamar [15]. Let I be the set of pixel sites (or voxels in three dimensions) in an image, and x be an image with one pixel value, xi at each site i 2 I . For a binary image,
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
xi 2 f0; 1g.
It is not possible to measure x directly, only a noisy approximation, y, where yi = xi + i , and the i are independent standard normal variables. At each site, we consider a neighborhood Ni of R 1 adjacent sites, including site i and centered at i. A pattern Æ is a particular coloring of the pixels in a neighborhood, formally, a vector of R zeros and ones. As the neighborhoods have R sites, there are 2R possible patterns. Bold face notation is used to represent a vector of variables in a neighborhood. Finally, Æ denotes the center pixel of a pattern. Let Q(Æ ) be the frequency of occurrence of any given pattern Æ in x. The probability mass function Q encapsulates some information about the local structure of the image. It can be shown that the estimate of xi , P Æ (y Æ)Q(Æ) x^i = PÆ i (1) ; Æ (yi Æ )Q(Æ ) where
(t) =
(2
2 f k2tk2 g
1 exp 2 )R=2
3
where xi is the vector of pixels in x in the neighborhood Ni , and w(xi ) is a weight used to emphasize the cost of the incorrect restoration of the relevant type of patterns. In order to derive the WMSE-optimal estimate, let the subscript J indicate the random selection of a site. Using such a random subscript, we may write W MSE = Ew(xJ )(xJ x^J )2 , and since for any random variables X , Y and Z and any function g we have
EZ X
E fZX jY g E fZ jY g
2
EZ (X g(Y ))2 ;
the WMSE-optimal estimate must be defined by x^i = g(yi ), where
g(t) =
E fw(xJ )xJ jyJ = tg : E fw(xJ )jyJ = tg
(4)
Following Meloche and Zamar [15], one easily obtains
P Æ (y Æ)w(Æ)Q(Æ) : x^i = PÆ i Æ (yi Æ )w(Æ )Q(Æ )
(5)
The WMSE criterion is clearly a generalization and yi is the vector of pixels in neighborhood Ni of of the AMSE criterion and the above derivation y, is the optimal average mean square error (AMSE) shows that the WMSE-optimal estimator is simply estimate. That is, it minimizes the objective function the AMSE-optimal estimator with Q(Æ ) replaced by X w(Æ)Q(Æ). Equation (4) makes no assumption about 1 2 AMSE = (2) the noise distribution and is general in this respect. jI j i2I E (xi x^i) ; Like Equation (1), Equation (5) assumes Gaussian over all estimates of the form xi = h(yi ) for some noise and varies for other distributions. The weighted probability w(Æ )Q(Æ ) is, in prach : RR ! R. tice, the same as the dictionary, D , in the relaxation Equation (1) is optimal in the sense of minimizing the AMSE. In the case where an object of interest method of Hancock and Kittler [9]. In our terminoloccupies only a small part of the image, however, the ogy, D = fÆ; w(Æ )Q(Æ ) 6= 0g. The above analysis AMSE criterion is a bad one. For example, for an au- shows that the dictionary arises as a natural consetomated inspection system that is designed to detect quence of the WMSE-optimal estimator. Selection the presence of a small object, a completely blank es- of neighborhood configurations for w(Æ )Q(Æ ) 6= 0 timate would result in a small AMSE, regardless of for thin objects is discussed in Section 3. Equations (1) and (5) give a simple algorithm to the image content. To make the restoration sensitive to rare but sig- compute an optimal restoration from y in a single nificant local phenomena, we introduce the weighted step. Improved results may be obtained by applying the restoration iteratively. We use the iterated conmean square error criterion given by ditional expectation (ICE) algorithm described by 1 X Meloche and Zamar [15] to improve restoration by 2 W MSE = jI j i2I Ew(xi )(xi x^i) ; (3) drawing support from sites further away.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
3 Algorithmic Considerations To create D for restoration of thin lines in twodimensional images we assume that neighborhoods are either all zeros or contain a single line segment. Neighborhoods containing a line segment all occur with the same weighted frequency, and all other neighborhoods have zero weight. When the size of a neighborhood is small, it is possible to enumerate manually the members of D . Table 1 summarizes D for a three-by-three neighborhood, R = 32 , and 29 = 512 possible neighborhoods. Note that D omits neighborhoods containing line segment terminations. and 1 are the weighted probabilities that a neighborhood is all blank or contains a line segment respectively. As the neighborhood size increases, the size of D precludes manual enumeration and an automated method becomes necessary. Our algorithm is described in [2]. While there are only 37 neighborhoods in the dictionary of Table 1 and the restoration algorithm is fast, as neighborhood size and the number of dimensions increase the dictionary grows rapidly and the restoration becomes too slow to be practical. If we precompute values of yi 0 and yi 1, then most of the computation lies in the summation of the norms kÆ yik2 which is repeated for every Æk 2 D. To make restoration of images with high noise levels practical requires that this process be as efficient as possible. One option is to capitalize on the fact that most of the Æk are nearly identical to another entry in D . This means that one can usually compute kÆk yi k2 quickly from kÆl yi k2 , l 6= k , Æl 2 D . To optimize this approach, consider the dictionary to be a completely connected graph where each Æk is a vertex and the cost of edge ek;l is the number of sites that differ between Æk and Æl , l 6= k . Given this representation of the dictionary, it is possible to compute a minimum-spanning tree using, for example, Prim’s algorithm [3]. Once the tree is constructed, it is then possible to compute all kÆk yi k2 in a minimum amount of time using a depth-first traversal. This approach leads to significant savings in computation by
4
avoiding redundant operations [2].
4 Results Figures 2 (a) and (b) show two-dimensional test images created by adding Gaussian noise to a binary image of thin objects (an arc of varying radius and a straight line). Application of ICE (three iterations) using the dictionary defined in Table 1, and = 0:5, yields the restored images in Figures 2 (c) and (d). Although the result in Figure 2(c) is trivial, a simple threshold would produce nearly the same result, the restoration in Figure 2(d) is also accurate and it is evident that the restoration is successful. Figures 2 (e) and (f) illustrate the effect of on the quality of the restoration. When is too low (Figure 2(e)), the algorithm is less likely to determine that a neighborhood contains all zeros and hallucinates lines where none exist. If is too high (Figure 2(f)), then the algorithm too easily assumes that a neighborhood is all zeros and misses the true presence of lines. One can exploit the contour and surface inference methods of Guy and Medioni [6, 7] for postprocessing by using a high value for . This gives a sparse set of points suitable for postprocessing, but avoids spurious background points. These results also illustrate the value of using w(Æ)Q(Æ) instead of Q(Æ). The sites in the test image that contain a line segment occupy a tiny portion of the image area, corresponding to > 0:95. Clearly, the trend depicted in Figure 2 shows that this value of would fail to restore the image. When critical parts of the image occur infrequently, the AMSE criterion is useless because it incurs only a small penalty by not restoring these parts. The WMSE criterion induces a proper restoration by increasing the cost of wrong restoration of these infrequent parts. Figure 3(a) shows a three-dimensional test pattern containing a single thin sheet of density one on a background of zeros. Points in the plot indicate voxels with a value greater than 0.5. The corresponding pattern with additive Gaussian noise, = 0:8, is shown in Figure 3(b). A single iteration of the binary restoration method using a neighborhood size of R = 53 , one iteration, and = 0:7 yields the
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
Æ
Class
5
Æ w(Æ)Q(Æ)
Blank
0
Through Center
1
1 36
Off Center
0
1 36
Table 1: Summary of D = fÆ; w(Æ )Q(Æ ) 6= 0g for a three-by-three neighborhood based on the assumption that neighborhoods are either blank or contain a line segment.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 2: Two-dimensional examples. (a) Test image ( = 0:1). (b) Test image ( = 0:4). (c) ICE restoration of (a) ( = 0:5). (d) ICE restoration of (b) ( = 0:5). (e) ICE restoration of (b) ( = 0:2). (f) ICE restoration of (b) ( = 0:8). Image size is 128 by 128. The restored images are thresholded so that pixels below 0.5 appear as zero and others as one.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98. restoration in Figure 3(c). The dictionary used for the restoration has 896 entries. A single iteration restores the sheet well, producing few errors. The images in Figure 3(d) and (e) give an impression of the noise levels involved. It is impossible to discern the thin sheet in Figure 3(d) but the cross section of the sheet in 3(e) is obvious. Figures 3(f) and (g) show what happens when we violate the noise model. The restoration in Figure 3(f) is the result of having correlated errors. We introduced correlations into the errors by smoothing the noise with a Gaussian kernel before adding it to the test pattern. The poor quality of the restoration occurs because correlated noise has structure that is mistaken for thin sheets. In Figure 3(g) the restoration is the result of having non-Gaussian errors. We contaminate the Gaussian noise by sampling from a gamma distribution for 10% of the voxels (randomly selected). The gamma distribution is scaled to have the same variance as the Gaussian, but its mean is non-zero. As a result, the errors are biased and asymmetric. We attribute the poor restoration primarily to this bias. In tomographic data, noise occurs as a result of limited photon statistics. Poisson noise at the photon sensors is transformed through the tomographic reconstruction process, resulting in an image with noise that is Gaussian, but varies with the voxel position and the contents of the image [11]. To apply the restoration method to tomographic data, we estimate the noise in the reconstructed image, and then assume it to be stationary. In general, tomographic images are not binary, materials in the world have widely varying densities. However, if the tomographic scanner is sensitive to some unique property of a material [16, 18], then it is still reasonable to assume a binary image; a material either has this property or it does not. Figures 3(h) and (i) show restorations of two synthetic tomographic images of a thin sheet. The synthetic data is derived from a simulation that emulates the scanner described by Sredniawski [19]. Both the photon statistics of the gamma source and the efficiency and spatial resolution of the detectors are modeled. The reconstruction computes ray
6
sums from the photon counts and then uses a filtered back projection method to produce images of volume slices. The thin sheet is based on the density properties of a typical plastic explosive, shaped into a thin sheet to confound detection. The thickness of the sheet and the scanner resolution are approximately equal. Figure 3(h) shows the successful restoration using a single iteration for an image with estimated noise of 0:5 and = 0:7. The restoration is still successful at 0:75 and = 0:9, Figure 3(i). Although a few errors are evident, the quality of restoration is sufficient for reliable detection.
5 Conclusions We have extrapolated the binary image restoration method of Meloche and Zamar to give a WMSE estimate of the data. The method is equivalent to weighting the true frequencies of possible neighborhoods. It successfully restores images with rare yet significant features. Our experimental results show the successful application of the method to both two- and threedimensional synthetic data, including tomographic reconstructions. The results also illustrate the importance of the parameter , the relative significance of background to object. At present we have no method to estimate what should be, so it requires manual tuning. The success of the method is encouraging and suggests that detection of thin concealed objects is possible even when noise levels are high and resolution is limited. Application to other difficult inspection tasks such as nondestructive testing for small stress fractures, and disbonds and delaminations in composite materials should be possible.
References [1] A. Blake and A. Zisserman. Visual Reconstruction. MIT Press, Cambridge, MA, 1987. [2] J. E. Boyd and J. Meloche. Binary restoration of thin objects in multidimensional imagery.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
Z
7
Z
Z
30
30
30
20
20
20
10
10
10
0
0
0
Y
Y
X
Y
X
(a)
X
(b)
(c)
Z
30
20
10
0 30
20 0 10
10 X
30
(d)
(e)
20
0
(f)
Z
Z
30
Y
20
Z
50
50
40
40
30
30
20
20
10
10
0
0
10
0 30
20 0 10
10 X
Y
Y
Y
20 30
(g)
0
X
X
(h)
(i)
Figure 3: Three-dimensional examples. (a) Test image containing a single thin sheet. (b) Image with additive Gaussian noise, = 0:8. (c) Restored image, = 0:8, = 0:7, 53 neighborhood. (d) Image slice in y z plane through unrestored image, = 0:8. (e) Restored image slice in y z plane, = 0:8, = 0:7, 53 neighborhood. (f) Restored image, = 0:8, noise is correlated by smoothing with a Gaussian kernel, smooth = 0:6. (g) Restored image = 0:8, mixture of Gaussian and gamma noise. (h) Restored tomographic image, 0:65, = 0:7, 53 neighborhood. (i) Restored tomographic image, 0:75, = 0:9, 53 neighborhood. Plot points indicate voxels with a value greater than 0.5.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
8
Technical Report 175, Department of Statistics, [12] K. Keeler. Map representations and codingbased priors for segmentation. In IEEE ConferUniversity of British Columbia, 1998. ence on Computer Vision and Pattern Recogni[3] G. Brassard and P. Bratley. Algorithmics Thetion 1991, pages 420–425, 1991. ory and Practice. Prentice Hall, Englewood [13] Y.G. Leclerc. Constructing simple stable deCliffs, NJ, 1988. scriptions for image partitioning. International [4] D. Geiger and F. Girosi. Parallel and determinJournal of Computer Vision, 3:73–102, 1989. istic algorithms from mrf’s: surface reconstruction. IEEE Transactions on Pattern Analysis [14] J.A. McDonald and A.B. Owen. Smoothing with split linear fits. Technometrics, 28(3):195– and Machine Intelligence, 13(5):401–412, May 208, August 1986. 1991. [15] J. Meloche and R.H. Zamar. Binary image [5] S. Geman and D. Geman. Stochastic relaxation, restoration. The Canadian Journal of Statistics, gibbs distributions, and the bayesian restoration 22(3):335–355, 1994. of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(2):721–741, [16] R.E. Morgado, B.Arnone, C.C. Cappiello, S.D. November 1984. Gardner, C.L. Hollas, L.E. Ussery, J.M. White, J.D. Zahrt, and R.A. Krauss. Prototype explo[6] G. Guy and G. Medioni. Inferring global persives detection system based on nuclear resoceptual contours from local features. Internance absorption in nitrogen. Proceedings of national Journal of Computer Vision, 20(1the SPIE, 2092:492–502, 1993. 2):113–133, October 1996. [17] P. Perona. Deformable kernels for early vi[7] G. Guy and G. Medioni. Inference of surfaces, sion. IEEE Transactions on Pattern Analysis 3d curves, and junctions from sparse noisy 3d and Machine Intelligence, 17(5):488–489, May data. IEEE Transactions on Pattern Analysis 1995. and Machine Intelligence, 19(11):1265–1277, [18] T. Seed, J.D Zahrt, and B.L. Berman. A virtual November 1997. prototype for an explosives detection system. In [8] P. Hall and D.M. Titterington. Edge-preserving Proceedings of the SPIE, volume 2092, pages and peak-preserving smoothing. Technomet482–491, 1993. rics, 34(4):429–440, November 1992. [19] J. Sredniawski. Detecting concealed explosives with gamma rays. The Industrial Physicist, [9] E.R. Hancock and J. Kittler. Edge-labeling uspages 24–27, March 1997. ing dictionary-based relaxation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(2):165–181, February 1990. [10] J.G. Harris, C. Koch, E. Staats, and J. Luo. Analog hardware for detecting discontinuities in early vision. International Journal of Computer Vision, 4:211–223, 1990. [11] A.C. Kak and M. Slaney. Principles of Computerized Tomographic Imaging. IEEE Press, New York, NY, 1988.
Abstract
gado et al. [16]. Therefore, detection of concealed explosives poses a particularly difficult problem, reWe present a method for restoration of noisy tomo- quiring the ability to identify thin objects in images graphic images for detecting thin objects, such as ex- with high noise and limited resolution. plosives. Use of a weighted mean-square estimate The goal of the work presented here is to find a optimizes the solution to place emphasis on the infresolution to this problem. Thin objects in two dimenquent, but significant local structure associated with sions are lines or curves no more than one pixel wide, thin objects. Experimental results show successful and in three dimensions they are sheets no more than restoration at very high noise levels. one voxel thick. Ultimately, we wish to detect these thin objects in data where the standard deviation of the noise is one to two times the density of the ob1 Introduction jects of interest. An important feature of the underlying image is its binary nature. Tomographic images Detection of certain objects in multidimensional improvide spatial density maps, so that, within a linear age data is a common problem in industrial inspectransformation, the density of an object is one on a tion. The objects may be a normal part of the subject background of zero. The assumption that there are of the inspection, or they may be anomalies. Deteconly two possible density values that can be mapped tion is more difficult if the object has been deliberonto ones and zeros is significant and, if exploited ately concealed, which is likely with explosives or as shown here, leads to successful detection in very illegal drugs. A common method used to conceal noisy data. While this work focuses on detection of malleable material is to form it into a thin sheet. It is concealed explosives, the method has potential for well known that this is effective for conventional Xapplication in other inspection tasks, such as stress ray imagery where a thin object is almost invisible fractures, and disbonds and delaminations in comunless viewed on end. The use of tomography ofposite materials. fers a solution to the problem of imaging thin sheets by viewing the object from many angles, but tomoWe view the detection problem as one of image graphic systems sensitive to the presence of explo- restoration. That is, we assume that if the noisy data sives are subject to high noise levels and limited res- are restored to a nearly noise free state, then the deolution. For example, see Seed et al. [18] and Mor- tection problem becomes trivial. Such is the case
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
(a)
(b)
Figure 1: Comparison of wide and thin object detection: (a) signal (solid) and smoothed signal (dashed) for a wide object, and (b) signal and smoothed signal for thin object.
when the imaging modality is sensitive to some characteristic unique to the objects of interest [16, 18]. The simplest method of restoration is to smooth the data with an isotropic linear shift-invariant lowpass filter. Such an approach tends to blur discontinuities, which is unacceptable for detection of thin objects. Figure 1 illustrates this fundamental problem. Smoothing a wide object blurs its edges but the object is still clearly present. The thin object all but disappears after smoothing. Therefore, isotropic smoothing of the data is not viable. Another option is to use many non-isotropic low-pass filters, each tuned to smooth only in one orientation. The required number of filters depends on the spatial support of the filters. This method requires much more computation and leaves the problem of reconciling the output of the different filters. Perona [17] describes steerable filters that are linear combinations of several non-isotropic filters. The number of filters required for accurate results increases with the aspect ratio of the individual isotropic filters. In the case of thin objects, the aspect ratio is high and a large number of filters is required. Several methods of restoration consider the presence of discontinuities in order to avoid smoothing across them. Geman and Geman [5], in their seminal work, model an image as a Markov random field (MRF) and include edge processes to mark discontinuities and avoid blurring. Blake and Zisserman [1] proposed the graduated non-convexity (GNC) algorithm that approximates an image with discontinuities, later implemented in hardware by Harris et al. [10]. Geiger and Girosi [4] derive GNC from MRF’s. Leclerc [13] optimizes a restoration by minimizing
2
the length of a code that describes the image regions. Keeler [12] proposes a similar method that incorporates topological considerations. These methods all rely, at some point, on a difference technique to infer the presence of a discontinuity or edge. As a result, these methods are effective where noise levels are not too great, but when the noise levels are high, reliable detection of edges is nearly impossible and these methods fail. Mcdonald and Owen [14] and Hall and Titterington [8] avoid reliance on edge detection by using split filters. The method presumes the existence of a specific number of discontinuities in each neighborhood of the data, and then computes the optimum positions of discontinuities to minimize the error between the data and the between-discontinuity fits. In multiple dimensions, the filter optimizes by performing an exhaustive search over a large number of possible discontinuity orientations and positions, and is very time-consuming. We present a method for binary image restoration, based on Meloche and Zamar [15], to exploit the binary nature of some data and achieve superior restoration results. For any neighborhood in a binary image, there is only a finite number of possible configurations of ones and zeros. Binary image restoration encapsulates knowledge of the structure of the true image in the function Q, the frequency distribution over the set of possible neighborhood configurations. In effect, selection of Q tunes the restoration to the characteristics of the objects to be found. Binary restoration is applicable to the detection of thin sheets so long as the underlying image is binary in nature. Our results indicate that the method is effective in restoring thin objects in noisy data, and that it is not necessary to know the true Q, but it is better to use a weighted Q, based on prior information.
2 Binary Image Restoration The following is a summary of the binary image restoration method proposed by Meloche and Zamar [15]. Let I be the set of pixel sites (or voxels in three dimensions) in an image, and x be an image with one pixel value, xi at each site i 2 I . For a binary image,
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
xi 2 f0; 1g.
It is not possible to measure x directly, only a noisy approximation, y, where yi = xi + i , and the i are independent standard normal variables. At each site, we consider a neighborhood Ni of R 1 adjacent sites, including site i and centered at i. A pattern Æ is a particular coloring of the pixels in a neighborhood, formally, a vector of R zeros and ones. As the neighborhoods have R sites, there are 2R possible patterns. Bold face notation is used to represent a vector of variables in a neighborhood. Finally, Æ denotes the center pixel of a pattern. Let Q(Æ ) be the frequency of occurrence of any given pattern Æ in x. The probability mass function Q encapsulates some information about the local structure of the image. It can be shown that the estimate of xi , P Æ (y Æ)Q(Æ) x^i = PÆ i (1) ; Æ (yi Æ )Q(Æ ) where
(t) =
(2
2 f k2tk2 g
1 exp 2 )R=2
3
where xi is the vector of pixels in x in the neighborhood Ni , and w(xi ) is a weight used to emphasize the cost of the incorrect restoration of the relevant type of patterns. In order to derive the WMSE-optimal estimate, let the subscript J indicate the random selection of a site. Using such a random subscript, we may write W MSE = Ew(xJ )(xJ x^J )2 , and since for any random variables X , Y and Z and any function g we have
EZ X
E fZX jY g E fZ jY g
2
EZ (X g(Y ))2 ;
the WMSE-optimal estimate must be defined by x^i = g(yi ), where
g(t) =
E fw(xJ )xJ jyJ = tg : E fw(xJ )jyJ = tg
(4)
Following Meloche and Zamar [15], one easily obtains
P Æ (y Æ)w(Æ)Q(Æ) : x^i = PÆ i Æ (yi Æ )w(Æ )Q(Æ )
(5)
The WMSE criterion is clearly a generalization and yi is the vector of pixels in neighborhood Ni of of the AMSE criterion and the above derivation y, is the optimal average mean square error (AMSE) shows that the WMSE-optimal estimator is simply estimate. That is, it minimizes the objective function the AMSE-optimal estimator with Q(Æ ) replaced by X w(Æ)Q(Æ). Equation (4) makes no assumption about 1 2 AMSE = (2) the noise distribution and is general in this respect. jI j i2I E (xi x^i) ; Like Equation (1), Equation (5) assumes Gaussian over all estimates of the form xi = h(yi ) for some noise and varies for other distributions. The weighted probability w(Æ )Q(Æ ) is, in prach : RR ! R. tice, the same as the dictionary, D , in the relaxation Equation (1) is optimal in the sense of minimizing the AMSE. In the case where an object of interest method of Hancock and Kittler [9]. In our terminoloccupies only a small part of the image, however, the ogy, D = fÆ; w(Æ )Q(Æ ) 6= 0g. The above analysis AMSE criterion is a bad one. For example, for an au- shows that the dictionary arises as a natural consetomated inspection system that is designed to detect quence of the WMSE-optimal estimator. Selection the presence of a small object, a completely blank es- of neighborhood configurations for w(Æ )Q(Æ ) 6= 0 timate would result in a small AMSE, regardless of for thin objects is discussed in Section 3. Equations (1) and (5) give a simple algorithm to the image content. To make the restoration sensitive to rare but sig- compute an optimal restoration from y in a single nificant local phenomena, we introduce the weighted step. Improved results may be obtained by applying the restoration iteratively. We use the iterated conmean square error criterion given by ditional expectation (ICE) algorithm described by 1 X Meloche and Zamar [15] to improve restoration by 2 W MSE = jI j i2I Ew(xi )(xi x^i) ; (3) drawing support from sites further away.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
3 Algorithmic Considerations To create D for restoration of thin lines in twodimensional images we assume that neighborhoods are either all zeros or contain a single line segment. Neighborhoods containing a line segment all occur with the same weighted frequency, and all other neighborhoods have zero weight. When the size of a neighborhood is small, it is possible to enumerate manually the members of D . Table 1 summarizes D for a three-by-three neighborhood, R = 32 , and 29 = 512 possible neighborhoods. Note that D omits neighborhoods containing line segment terminations. and 1 are the weighted probabilities that a neighborhood is all blank or contains a line segment respectively. As the neighborhood size increases, the size of D precludes manual enumeration and an automated method becomes necessary. Our algorithm is described in [2]. While there are only 37 neighborhoods in the dictionary of Table 1 and the restoration algorithm is fast, as neighborhood size and the number of dimensions increase the dictionary grows rapidly and the restoration becomes too slow to be practical. If we precompute values of yi 0 and yi 1, then most of the computation lies in the summation of the norms kÆ yik2 which is repeated for every Æk 2 D. To make restoration of images with high noise levels practical requires that this process be as efficient as possible. One option is to capitalize on the fact that most of the Æk are nearly identical to another entry in D . This means that one can usually compute kÆk yi k2 quickly from kÆl yi k2 , l 6= k , Æl 2 D . To optimize this approach, consider the dictionary to be a completely connected graph where each Æk is a vertex and the cost of edge ek;l is the number of sites that differ between Æk and Æl , l 6= k . Given this representation of the dictionary, it is possible to compute a minimum-spanning tree using, for example, Prim’s algorithm [3]. Once the tree is constructed, it is then possible to compute all kÆk yi k2 in a minimum amount of time using a depth-first traversal. This approach leads to significant savings in computation by
4
avoiding redundant operations [2].
4 Results Figures 2 (a) and (b) show two-dimensional test images created by adding Gaussian noise to a binary image of thin objects (an arc of varying radius and a straight line). Application of ICE (three iterations) using the dictionary defined in Table 1, and = 0:5, yields the restored images in Figures 2 (c) and (d). Although the result in Figure 2(c) is trivial, a simple threshold would produce nearly the same result, the restoration in Figure 2(d) is also accurate and it is evident that the restoration is successful. Figures 2 (e) and (f) illustrate the effect of on the quality of the restoration. When is too low (Figure 2(e)), the algorithm is less likely to determine that a neighborhood contains all zeros and hallucinates lines where none exist. If is too high (Figure 2(f)), then the algorithm too easily assumes that a neighborhood is all zeros and misses the true presence of lines. One can exploit the contour and surface inference methods of Guy and Medioni [6, 7] for postprocessing by using a high value for . This gives a sparse set of points suitable for postprocessing, but avoids spurious background points. These results also illustrate the value of using w(Æ)Q(Æ) instead of Q(Æ). The sites in the test image that contain a line segment occupy a tiny portion of the image area, corresponding to > 0:95. Clearly, the trend depicted in Figure 2 shows that this value of would fail to restore the image. When critical parts of the image occur infrequently, the AMSE criterion is useless because it incurs only a small penalty by not restoring these parts. The WMSE criterion induces a proper restoration by increasing the cost of wrong restoration of these infrequent parts. Figure 3(a) shows a three-dimensional test pattern containing a single thin sheet of density one on a background of zeros. Points in the plot indicate voxels with a value greater than 0.5. The corresponding pattern with additive Gaussian noise, = 0:8, is shown in Figure 3(b). A single iteration of the binary restoration method using a neighborhood size of R = 53 , one iteration, and = 0:7 yields the
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
Æ
Class
5
Æ w(Æ)Q(Æ)
Blank
0
Through Center
1
1 36
Off Center
0
1 36
Table 1: Summary of D = fÆ; w(Æ )Q(Æ ) 6= 0g for a three-by-three neighborhood based on the assumption that neighborhoods are either blank or contain a line segment.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 2: Two-dimensional examples. (a) Test image ( = 0:1). (b) Test image ( = 0:4). (c) ICE restoration of (a) ( = 0:5). (d) ICE restoration of (b) ( = 0:5). (e) ICE restoration of (b) ( = 0:2). (f) ICE restoration of (b) ( = 0:8). Image size is 128 by 128. The restored images are thresholded so that pixels below 0.5 appear as zero and others as one.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98. restoration in Figure 3(c). The dictionary used for the restoration has 896 entries. A single iteration restores the sheet well, producing few errors. The images in Figure 3(d) and (e) give an impression of the noise levels involved. It is impossible to discern the thin sheet in Figure 3(d) but the cross section of the sheet in 3(e) is obvious. Figures 3(f) and (g) show what happens when we violate the noise model. The restoration in Figure 3(f) is the result of having correlated errors. We introduced correlations into the errors by smoothing the noise with a Gaussian kernel before adding it to the test pattern. The poor quality of the restoration occurs because correlated noise has structure that is mistaken for thin sheets. In Figure 3(g) the restoration is the result of having non-Gaussian errors. We contaminate the Gaussian noise by sampling from a gamma distribution for 10% of the voxels (randomly selected). The gamma distribution is scaled to have the same variance as the Gaussian, but its mean is non-zero. As a result, the errors are biased and asymmetric. We attribute the poor restoration primarily to this bias. In tomographic data, noise occurs as a result of limited photon statistics. Poisson noise at the photon sensors is transformed through the tomographic reconstruction process, resulting in an image with noise that is Gaussian, but varies with the voxel position and the contents of the image [11]. To apply the restoration method to tomographic data, we estimate the noise in the reconstructed image, and then assume it to be stationary. In general, tomographic images are not binary, materials in the world have widely varying densities. However, if the tomographic scanner is sensitive to some unique property of a material [16, 18], then it is still reasonable to assume a binary image; a material either has this property or it does not. Figures 3(h) and (i) show restorations of two synthetic tomographic images of a thin sheet. The synthetic data is derived from a simulation that emulates the scanner described by Sredniawski [19]. Both the photon statistics of the gamma source and the efficiency and spatial resolution of the detectors are modeled. The reconstruction computes ray
6
sums from the photon counts and then uses a filtered back projection method to produce images of volume slices. The thin sheet is based on the density properties of a typical plastic explosive, shaped into a thin sheet to confound detection. The thickness of the sheet and the scanner resolution are approximately equal. Figure 3(h) shows the successful restoration using a single iteration for an image with estimated noise of 0:5 and = 0:7. The restoration is still successful at 0:75 and = 0:9, Figure 3(i). Although a few errors are evident, the quality of restoration is sufficient for reliable detection.
5 Conclusions We have extrapolated the binary image restoration method of Meloche and Zamar to give a WMSE estimate of the data. The method is equivalent to weighting the true frequencies of possible neighborhoods. It successfully restores images with rare yet significant features. Our experimental results show the successful application of the method to both two- and threedimensional synthetic data, including tomographic reconstructions. The results also illustrate the importance of the parameter , the relative significance of background to object. At present we have no method to estimate what should be, so it requires manual tuning. The success of the method is encouraging and suggests that detection of thin concealed objects is possible even when noise levels are high and resolution is limited. Application to other difficult inspection tasks such as nondestructive testing for small stress fractures, and disbonds and delaminations in composite materials should be possible.
References [1] A. Blake and A. Zisserman. Visual Reconstruction. MIT Press, Cambridge, MA, 1987. [2] J. E. Boyd and J. Meloche. Binary restoration of thin objects in multidimensional imagery.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
Z
7
Z
Z
30
30
30
20
20
20
10
10
10
0
0
0
Y
Y
X
Y
X
(a)
X
(b)
(c)
Z
30
20
10
0 30
20 0 10
10 X
30
(d)
(e)
20
0
(f)
Z
Z
30
Y
20
Z
50
50
40
40
30
30
20
20
10
10
0
0
10
0 30
20 0 10
10 X
Y
Y
Y
20 30
(g)
0
X
X
(h)
(i)
Figure 3: Three-dimensional examples. (a) Test image containing a single thin sheet. (b) Image with additive Gaussian noise, = 0:8. (c) Restored image, = 0:8, = 0:7, 53 neighborhood. (d) Image slice in y z plane through unrestored image, = 0:8. (e) Restored image slice in y z plane, = 0:8, = 0:7, 53 neighborhood. (f) Restored image, = 0:8, noise is correlated by smoothing with a Gaussian kernel, smooth = 0:6. (g) Restored image = 0:8, mixture of Gaussian and gamma noise. (h) Restored tomographic image, 0:65, = 0:7, 53 neighborhood. (i) Restored tomographic image, 0:75, = 0:9, 53 neighborhood. Plot points indicate voxels with a value greater than 0.5.
Boyd and Meloche, PAMI, Vol. 20, No. 6, p647–651, June 98.
8
Technical Report 175, Department of Statistics, [12] K. Keeler. Map representations and codingbased priors for segmentation. In IEEE ConferUniversity of British Columbia, 1998. ence on Computer Vision and Pattern Recogni[3] G. Brassard and P. Bratley. Algorithmics Thetion 1991, pages 420–425, 1991. ory and Practice. Prentice Hall, Englewood [13] Y.G. Leclerc. Constructing simple stable deCliffs, NJ, 1988. scriptions for image partitioning. International [4] D. Geiger and F. Girosi. Parallel and determinJournal of Computer Vision, 3:73–102, 1989. istic algorithms from mrf’s: surface reconstruction. IEEE Transactions on Pattern Analysis [14] J.A. McDonald and A.B. Owen. Smoothing with split linear fits. Technometrics, 28(3):195– and Machine Intelligence, 13(5):401–412, May 208, August 1986. 1991. [15] J. Meloche and R.H. Zamar. Binary image [5] S. Geman and D. Geman. Stochastic relaxation, restoration. The Canadian Journal of Statistics, gibbs distributions, and the bayesian restoration 22(3):335–355, 1994. of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(2):721–741, [16] R.E. Morgado, B.Arnone, C.C. Cappiello, S.D. November 1984. Gardner, C.L. Hollas, L.E. Ussery, J.M. White, J.D. Zahrt, and R.A. Krauss. Prototype explo[6] G. Guy and G. Medioni. Inferring global persives detection system based on nuclear resoceptual contours from local features. Internance absorption in nitrogen. Proceedings of national Journal of Computer Vision, 20(1the SPIE, 2092:492–502, 1993. 2):113–133, October 1996. [17] P. Perona. Deformable kernels for early vi[7] G. Guy and G. Medioni. Inference of surfaces, sion. IEEE Transactions on Pattern Analysis 3d curves, and junctions from sparse noisy 3d and Machine Intelligence, 17(5):488–489, May data. IEEE Transactions on Pattern Analysis 1995. and Machine Intelligence, 19(11):1265–1277, [18] T. Seed, J.D Zahrt, and B.L. Berman. A virtual November 1997. prototype for an explosives detection system. In [8] P. Hall and D.M. Titterington. Edge-preserving Proceedings of the SPIE, volume 2092, pages and peak-preserving smoothing. Technomet482–491, 1993. rics, 34(4):429–440, November 1992. [19] J. Sredniawski. Detecting concealed explosives with gamma rays. The Industrial Physicist, [9] E.R. Hancock and J. Kittler. Edge-labeling uspages 24–27, March 1997. ing dictionary-based relaxation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(2):165–181, February 1990. [10] J.G. Harris, C. Koch, E. Staats, and J. Luo. Analog hardware for detecting discontinuities in early vision. International Journal of Computer Vision, 4:211–223, 1990. [11] A.C. Kak and M. Slaney. Principles of Computerized Tomographic Imaging. IEEE Press, New York, NY, 1988.
