Project 2: Where are My Glasses?
The image above is a blurred version of a picture that I took at the local SPCA. As we will see, a naive approach to de-blurring produces garbage, a characteristic feature of ill-posed problems. In order to reconstruct the image, we need to regularize the reconstruction, just as we would regularize an ill-posed least squares problem. We will use Tikhonov regularization, as described in class. However, this involves choosing the value of a regularization parameter. Your mission is to investigate the dependence on the parameter, and to investigate three approaches to choosing this parameter (mostly) automatically.
You are given the file blurry.png (the blurred image shown above). You are responsible for the following deliverables, which can either by submitted as a Jupyter notebook (recommended) or in standalone code:
Codes to compute "optimal" values of the regularization parameter $\lambda$ via the discrepancy principle, the L-curve, and generalized cross-validation.
A report that addresses the questions posed in the rest of this prompt.
This project is inspired by the project on image deblurring by James G. Nagy and Dianne P. O'Leary (Image Deblurring: I Can See Clearly Now) in Computing in Science and Engineering; Project: Vol. 5, No. 3, May/June 2003, pp. 82-85; Solution: Vol. 5, No. 4, July/August 2003). Other useful references include:
Hansen, Nagy, O'Leary. Deblurring Images: Matrices, Spectra, and Filtering, SIAM 2006.
Hansen and O'Leary. "The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems" in SIAM J. Sci. Comput., Vol 14, No. 6, November 1993.
Golub, Heath, and Wahba. "Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter" in Technometrics, Vol. 21, No. 2, May 1979.
Golub and von Matt. "Generalized Cross-Validation for Large-Scale Problems" in Journal of Computational and Graphical Statistics, Vol. 6, No. 1, March 1997.
Morozov. Methods for Solving Incorrectly Posed Problems, Springer-Verlag, 1984.
It should be possible to do this assignment based only on ideas presented in this prompt and in the lectures. However, you are always welcome to use any ideas or code you find in the literature, including in the references noted above, provided you give an appropriate citation.
Julia setup
We use the following Julia packages for this project.
FileIO: The FileIO library provides extensible support for loading different types of files; we use it to work with PNG files.
Images: The Images package provides various tools for working with images.
FFTW: The FFTW package provides Fourier transform tools.
Colors: Manipulates different color spaces
Roots: Solve one-dimensional nonlinear equations
Plots: Create plots in Julia
You will also want to install the ImageMagick or QuartzImageIO (on Mac) to work with FileIO to load images. You do not need to include these packages directly with using, but they will be loaded by the FileIO library.
Once we have the prerequisites, we can load the blurred image. If you have the data locally, it is fine to just use load("blurry.png") to get the image; but one can also combine the load with a download command to fetch the data from an online source.
Image array manipulation
Images are represented in Julia as two-dimensional arrays of color values. In the default representation of the PNG file that we downloaded, each pixel color is represented by 24 bits, with 8 bits each for the red, green, and blue intensity values (RGB). Each 8 bit number is represented as a fixed-point number with values ranging from 0.0 = 0/255 to 1.0 = 255/255 (a Normed{Uint8,8} type in Julia).
Matrix{RGB{N0f8}} (alias for Array{RGB{Normed{UInt8, 8}}, 2})The picture was generated by a blurring operation that works on each color plane independently. To deblur, we want to work on the arrays for one color plane at a time. Therefore, we will manipulate the data in terms of three floating-point arrays (one each for the red, green, and blue levels), which we construct via the channelview function. The colorview function maps these floating point arrays back to a viewable image. We define a convenience function arrays_to_img that does this conversion and makes sure that all the entries of the arrays are clipped back to the interval from 0 to 1.
img_rgb_pixels (generic function with 1 method)arrays_to_img (generic function with 1 method)Now we can compute the 3-array representation of the image and illustrate the helpers.
3×249×320 reinterpret(reshape, Float32, ::Array{RGB{Float32},2}) with eltype Float32:
[:, :, 1] =
0.360784 0.345098 0.34902 0.352941 … 0.376471 0.388235 0.396078 0.388235
0.443137 0.435294 0.443137 0.447059 0.447059 0.458824 0.470588 0.466667
0.545098 0.541176 0.54902 0.556863 0.533333 0.54902 0.560784 0.560784
[:, :, 2] =
0.360784 0.345098 0.34902 0.352941 … 0.392157 0.403922 0.407843 0.392157
0.439216 0.431373 0.439216 0.447059 0.458824 0.470588 0.478431 0.466667
0.541176 0.537255 0.54902 0.556863 0.541176 0.552941 0.564706 0.560784
[:, :, 3] =
0.356863 0.341176 0.34902 0.352941 … 0.407843 0.415686 0.415686 0.392157
0.435294 0.431373 0.439216 0.447059 0.470588 0.482353 0.482353 0.462745
0.537255 0.533333 0.54902 0.556863 0.545098 0.560784 0.568627 0.556863
;;; …
[:, :, 318] =
0.364706 0.352941 0.352941 0.352941 … 0.329412 0.341176 0.352941 0.364706
0.447059 0.443137 0.447059 0.45098 0.407843 0.423529 0.435294 0.447059
0.552941 0.552941 0.556863 0.560784 0.509804 0.52549 0.541176 0.552941
[:, :, 319] =
0.364706 0.352941 0.34902 0.352941 … 0.345098 0.356863 0.368627 0.376471
0.447059 0.439216 0.447059 0.45098 0.419608 0.435294 0.447059 0.454902
0.552941 0.54902 0.552941 0.556863 0.517647 0.533333 0.54902 0.556863
[:, :, 320] =
0.364706 0.34902 0.34902 0.352941 … 0.360784 0.372549 0.384314 0.384314
0.447059 0.439216 0.443137 0.447059 0.435294 0.447059 0.458824 0.462745
0.54902 0.545098 0.552941 0.556863 0.52549 0.541176 0.556863 0.560784npixels = 79680