Discontinuous Regression and Image Processing
Jump Regression Models
Observed digital images often contain noise which could be incurred from different sources. For example, noise can come from the digitization process described above. Noise can also be produced in poor lighting situations where images are taken. Assuming that the noise is additive to its true intensity values, an image can be described by the following two-dimensional regression model: Zi, j = f(xi, yj) + εi, j, (xi, yj) ∈ Ω, i = 1, …, n1, j = 1, …, n2, where f is the true image intensity function, εi, j denotes the noise at the (i, j)-th pixel, Zi, j is the observed intensity value at the (i, j)-th pixel, and Ω is the design space. The image intensity function f often has discontinuities. For instance, the boundary curves of objects in a photograph are positions where f has jumps. Because our human-eye systems have evolved to make use of the boundary curves for recognizing objects, jumps in f are important features of the image. In image processing, jump locations in f are called step edges, and jump locations in the first-order derivatives of f are called roof edges. Therefore, edge detection and edge-preserving image restoration in image processing are essentially the same tasks as jump detection and jump-preserving surface estimation in jump regression.
Package Overview
The DRIP package provides functionality for three categories of image analysis problems: edge detection, edge-preserving noise removal and blind image deblurring. Edge detection is often used in image data compression. Noise removal and image deblurring are important for improving human and machine perception of the images. The package is designed for general image analysis without restrictions on the specific image type. To use DRIP in practice, the workflow often involves parameter tuning before model estimation.
Step Edge Detection
The step edge detector in DRIP is implemented by the
stepEdge() function. We apply it to the SAR image which is
included in the package. The output of stepEdge() is a
matrix of 0’s and 1’s. It is worth noting that the bandwidth is
specified in terms of the number of pixels. Also, we have exchanged 0’s
and 1’s in the visualization with the image() function in
order to have the jump points shown in black.
stepedge <- stepEdge(image = sar, bandwidth = 10, thresh = 17,
degree = 1)
par(mfrow = c(1, 2), mar = c(3, 1, 1, 1), xaxt = "n", yaxt = "n")
image(sar, col = gray(c(0:255)/255))
image(1 - stepedge, col = gray(c(0:255)/255))
Two parameters, bandwidth and threshold, need to be specified by the
user. The stepEdgeParSel() function implements a
bootstrap-based data-driven procedure that helps determine these two
parameter values. In the following illustration with the SAR image, we
provide two bandwidths and two thresholds and select the best
combination using stepEdgeParSel(). Since this is a
bootstrap procedure, the number of bootstrap samples needs to be given,
and a random seed is required to ensure reproducibility. The function
returns a matrix of d̂KQ, an
edge detection performance measure, for each combination. It also
reports the selected bandwidth and threshold parameters associated with
the smallest d̂KQ
value.
set.seed(24)
parSel <- stepEdgeParSel(image = sar, bandwidth = c(9, 10), degree = 1,
thresh = c(17, 21), nboot = 10)
print(parSel, type = "all")
#> The bootstrap matrix:
#> thresh=17 thresh=21
#> bandwidth=9 0.009777572 0.012336507
#> bandwidth=10 0.008982834 0.011186066
#> The selected bandwidth: 10
#> The selected threshold: 17Since the bandwidth specifies the neighborhood size for local
smoothing, positive integers in the range of [1, 20] are usually reasonable candidates. It
is more difficult to obtain a preliminary range for reasonable threshold
values. To handle this issue, one may use the stepDiff()
function, which computes the edge detection criterion at each pixel.
Since the number of edge pixels is small relative to the image
resolution, upper percentiles provide good clues about reasonable values
for the threshold parameter.
Image Denoising
The threeStage() function in the package implements a
three-stage approach for image denoising. The detected step edges are
supplied using the argument edge1. Notably, users can also
input the detected roof edges using the argument edge2. The
resulting estimator preserves both step and roof edges of the image. In
cases where roof edge detection is deemed unnecessary, users can simply
provide a matrix of zeros to edge2. The code below shows
the estimated SAR image after applying threeStage(). It can
be seen that the estimation preserves the discontinuities at places
where edges have been successfully detected. At places where the edge
detector fails to flag the edge points, the estimation still blurs the
jumps.
fit <- threeStage(image = sar, bandwidth = 4, edge1 = stepedge,
edge2 = array(0, dim(sar)))
par(mfrow = c(1, 1), mar = c(3, 1, 1, 1), xaxt = "n", yaxt = "n")
image(fit, col = gray(c(0:255)/255))
This three-stage approach requires a bandwidth parameter. It can be
selected by minimizing the leave-one-out cross validation score. The
threeStageParSel() function implements the cross validation
procedure.
Blind Image Deblurring
In addition to having noise, images may also have blur involved. For instance, in astronomical imaging, ground-based imaging systems are subject to blurring due to the rapidly changing index of refraction of the atmosphere. In photography, out-of-focus or camera shake often results in blurred images. In medical imaging, blurred x-rays or mammograms are almost inevitable because the medical imaging systems limit the intensity of the incident radiation in order to protect the patient’s health. In the image processing literature, a commonly used model to describe the relationship between the original image and its observed but blurred version is Z(xi, yj) = G{f}(xi, yj) + εi, j, where G{f} is the convolution between g and f defined by G{f}(x, y) = g ⊗ f(x, y) = ∫∫Ωg(x − u, y − v; x, y)f(u, v) dudv. Here g is called the point spread function. It describes how the original image is spatially degraded (i.e., blurred). In most references, it is further assumed that the blurring is location (or spatially) invariant. That is, g(u, v; x, y) does not depend on (x, y). Blind image deblurring is for estimating f from Z when the point spread function g is not completely specified.
The jpex() function implements a blind deblurring method
using the jump-preserving extrapolation (JPEX) technique. Here we apply
it to a blurry stop-sign image, which is shown below. The function
returns a list of two: deblurred, the deblurred image, and
edge, the detected blurry pixels. It can be seen that the
JPEX method is able to deblur the image well. The arguments
alpha and sigma specify the significance level
and noise level, respectively.
deblur <- jpex(image = stopsign, bandwidth = 2, sigma = 0.00623,
alpha = 0.001)
names(deblur)
#> [1] "deblurred" "edge"
par(mfrow = c(1, 2), mar = c(3, 1, 1, 1), xaxt = "n", yaxt = "n")
image(stopsign, col = gray(c(0:255)/255))
image(deblur$deblurred, col = gray(c(0:255)/255))
The jpex() function requires a bandwidth value and the
noise level. Both can be determined by the cv.jpex()
function.