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: \[
Z_{i, j} = f(x_i, y_j) + \varepsilon_{i, j}, \; (x_i, y_j) \in \Omega,
\; i = 1, \ldots, \; n_1, \; j = 1, \ldots, n_2,
\] where \(f\) is the true image
intensity function, \(\varepsilon_{i,
j}\) denotes the noise at the \((i,
j)\)-th pixel, \(Z_{i,j}\) is
the observed intensity value at the \((i,
j)\)-th pixel, and \(\Omega\) 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. DRIP has two functions for this purpose. The
first is edgeParSelPilot(). Users can specify several
preliminary bandwidth values. Since the bandwidth specifies the
neighborhood size for local smoothing, positive integers in the range of
\([1, 20]\) are usually reasonable
candidates. The edgeParSelPilot() function would plot the
image of the edge detection statistics for each given bandwidth. Users
can narrow down the choices based on the visual impression. In addition,
the function outputs a few upper quantiles of the edge detection
statistics for each bandwidth. Those quantiles offer clues to what
reasonable threshold values can be, as edges are usually sparse in an
image. In the following illustration with the SAR image, we provide
three bandwidths and three probabilities to
edgeParSelPilot().
edgeParSelPilot(sar, edgeType = "step", degree = 1, bandwidth = c(6, 8, 10), probs = c(0.75, 0.85, 0.95))
#> probs=0.75 probs=0.85 probs=0.95
#> bandwidth=6 10.688440 13.79875 20.18085
#> bandwidth=8 9.610952 12.59395 19.38523
#> bandwidth=10 9.226387 12.21794 19.25910
Bandwidth values around 10 seem reasonable as the visualization is
less noisy. The quantiles suggest that the threshold should be around
20. Next, users can use the second function,
stepEdgeParSel(), to finalize the parameter choice. The
stepEdgeParSel() function implements a bootstrap-based
data-driven procedure that estimates the edge detection performance for
each combination of the bandwidth and threshold. 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 \(\widehat{d}_{KQ}\), an edge detection
performance measure. It also reports the selected bandwidth and
threshold parameters associated with the smallest \(\widehat{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: 17
Image Denoising
The restore3Stage() 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 restore3Stage(). 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 <- restore3Stage(image = sar, bandwidth = 4, step_edge = stepedge,
roof_edge = 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
restore3StageParSel() function implements the cross
validation procedure.
bw_3stage <- restore3StageParSel(image = sar, bandwidth = 4:5,
step_edge = stepedge,
roof_edge = array(0, dim(sar)))
print(bw_3stage, type = "all")
#> The cross validation scores:
#> bandwidth=4 bandwidth=5
#> CV-score 293.0294 300.4995
#> The selected bandwidth: 4
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(x_i, y_j) = G\{f\}(x_i, y_j) + \varepsilon_{i, j},
\] where \(G\{f\}\) is the
convolution between \(g\) and \(f\) defined by \[
G\{f\}(x,y) = g\otimes f(x, y) = \int\int_\Omega 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.
cv.out <- cv.jpex(image = stopsign, bandwidths = c(2, 3), ncpus = 1)
print(cv.out, type = "all")
#> The cross validation scores:
#> [1] 5.591376e-05 1.421656e-04
#> The selected bandwidth: 2
#> The estimated sigma: 0.006231292