+
RX Anomaly Detector
+
The RX anomaly detector uses the squared Mahalanobis distance as a measure of
+how anomalous a pixel is with respect to an assumed background. The SPy
+rx function computes RX scores for an
+array of image pixels. The squared Mahalanobis distance is given by
+
+
+
where 
is the pixel spectrum, 
is the background
+mean, and 
is the background covariance [Reed_Yu_1990].
+
If no background statistics are passed to the rx function,
+background statistics will be estimated from the array of pixels for which the
+RX scores are to be calculated.
+
In [44]: rxvals = rx(img)
+
+
+
To declare pixels as anomalous, we need to specify a threshold RX score. For
+example, we could choose all image pixels whose RX score has a probability of
+less than 0.001 with respect to the background:
+
In [45]: from scipy.stats import chi2
+
+In [46]: nbands = img.shape[-1]
+
+In [47]: P = chi2.ppf(0.999, nbands)
+
+In [48]: v = imshow(1 * (rxvals > P))
+
+
+
+
Rather than specifying a threshold for anomalous pixels, one can also simply
+view an image of raw RX scores, where brighter pixels are considered “more anomalous”:
+
In [49]: v = imshow(rxvals)
+
+
+
+
For the sample image, only a few pixels are visible in the image of RX scores
+because a linear color scale is used and there is a very small number of pixels with RX
+scores much higher than other pixels. This is apparent from viewing the histogram
+of the RX scores.
+
In [50]: import matplotlib.pyplot as plt
+
+In [51]: f = plt.figure()
+
+In [52]: h = plt.hist(rxvals.ravel(), 200, log=True)
+
+In [53]: h = plt.grid()
+
+
+
+
The outliers are not obvious in the histogram, so let’s print their values:
+
In [54]: print(np.sort(rxvals.ravel())[-10:])
+[ 665.17211462 675.85801536 683.58190673 731.4872873 739.95211335
+ 906.98669373 956.49972325 1703.20957949 2336.11246149 9018.65517253]
+
+
+
To view greater detail in the RX image, we can adjust the lower and upper limits
+of the image display. Since we are primarily interested in the most anomalous
+pixels, we will set the black level to the 99th percentile of the RX values’
+cumulative histogram and set the white point to the 99.99th percentile:
+
In [55]: v = imshow(rxvals, stretch=(0.99, 0.9999))
+
+
+
+
We can see the new RGB data limits by inspecting the returned
+ImageView object:
+
In [56]: print(v)
+ImageView object:
+ Display bands : [0]
+ Interpolation : <default>
+ RGB data limits :
+ R: [291.6938067178437, 956.4997232540622]
+ G: [291.6938067178437, 956.4997232540622]
+ B: [291.6938067178437, 956.4997232540622]
+
+
+
Note that we could also have set the contrast stretch to explicit RX values
+(vice percentile values) by specifying the bounds keyword instead of stretch.
+
If an image contains regions with different background materials, then the
+assumption of a single mean/covariance for background pixels can reduce
+performance of the RX anomaly detector. In such situations, better results
+can be obtained by dynamically computing background statistics in a neighborhood
+around each pixel being evaluated.
+
To compute local background statistics for
+each pixel, the rx function accepts an optional window argument, which
+specifies an inner/outer window within which to calculate background statistics
+for each pixel being evaluated. The outer window is the window within which
+background statistics will be calculated. The inner window is a smaller window
+(within the outer window) indicating an exclusion zone within which pixels are
+to be ignored. The purpose of the inner window is to prevent potential anomaly/target
+pixels from “polluting” background statistics.
+
For example, to compute RX scores with background statistics computed from a
+21 
21 pixel window about the pixel being evaluated, with an exclusion window of
+5 5 pixels, the function would be called as follows:
+
In [57]: rxvals = rx(img, window=(5,21))
+
+
+
While the use of a windowed background will often improve results for images containing
+multiple background materials, it does so at the expense of introducing two
+issues. First, the sizes of the inner and outer windows must be specified such
+that the resulting covariance has full rank. That is, if 
and
+
represent the pixel widths of the inner and outer windows, respectively,
+then 
must be at least as large as the number of
+image bands. Second, recomputing the estimated background covariance for each
+pixel in the image makes the computational complexity of of the RX score
+computation orders of magnitue greater.
+
As a compromise between a fixed background and recomputation of mean & covariance
+for each pixel, rx can be passed a global covariance estimate in addition
+to the window argument. In this case, only the background mean within the window will
+be recomputed for each pixel. This significanly reduces computation time
+for the windowed algorithm and removes the size limitation on the window (except
+that the outer window must be larger than the inner). For example, since our sample image has ground
+cover classes labeled, we can compute the average covariance over those ground
+cover classes and use the result as an estimate of the global covariance.
+
In [58]: C = cov_avg(img, gt)
+
+In [59]: rxvals = rx(img, window=(5,21), cov=C)
+
+
+
+
Matched Filter
+
The matched filter is a linear detector given by the formula
+
+
+
where 
is the target mean, is the background
+mean, and is the covariance. The matched filter response is
+scaled such that the response is zero when the input is equal to the background
+mean and equal to one when the pixel is equal to the target mean. Like the
+rx function the SPy
+matched_filter function will estimate
+background statistics from the input image if no background statistics are
+specified.
+
Let’s select the image pixel at (row, col) = (8, 88) as our target, use a
+global background statistics estimate, and plot all pixels whose matched filter
+scores are greater than 0.2.
+
In [60]: t = img[8, 88]
+
+In [61]: mfscores = matched_filter(img, t)
+
+In [62]: v = imshow(1 * (mfscores > 0.2))
+
+
+
+
As with the rx function, matched_filter can be applied using
+windowed background statistics (optionally with a global covariance estimate).
+
+
See also
+
ace:
+
+Adaptive Coherence/Cosine Estimator (ACE)
+
+
+