Edit or run this notebook

Project 2: Where are My Glasses?

Due: 2022-03-25

blurred image

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 $λ$ 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.

md"""

# Project 2: Where are My Glasses?

## Due: 2022-03-25

![blurred image](https://github.com/dbindel/cs4220-s20/raw/master/hw/code/proj1/data/blurry.png)

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](http://ieeexplore.ieee.org/document/1196312/))

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*](http://epubs.siam.org/doi/book/10.1137/1.9780898718874), SIAM 2006.

- Hansen and O'Leary. ["The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems"](http://epubs.siam.org/doi/abs/10.1137/0914086) in *SIAM J. Sci. Comput.*, Vol 14, No. 6, November 1993.

35.6 ms

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.

md"""

## Julia setup

We use the following Julia packages for this project.

- [FileIO](https://github.com/JuliaIO/FileIO.jl): The FileIO library provides extensible support for loading different types of files; we use it to work with PNG files.

- [Images](https://juliaimages.org/latest/): The Images package provides various tools for working with images.

- [FFTW](https://github.com/JuliaMath/FFTW.jl): The FFTW package provides Fourier transform tools.

- [Colors](https://github.com/JuliaGraphics/Colors.jl): Manipulates different color spaces

- [Roots](https://github.com/JuliaMath/Roots.jl): Solve one-dimensional nonlinear equations

- [Plots](http://docs.juliaplots.org/latest/): Create plots in Julia

You will also want to install the [ImageMagick](https://github.com/JuliaIO/ImageMagick.jl) or [QuartzImageIO](https://github.com/JuliaIO/QuartzImageIO.jl) (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.

"""

713 μs

using Images

2.7 s

using FileIO

161 μs

using FFTW

150 μs

using Colors

171 μs

using Roots

62.2 ms

using Plots

2.5 s

using DelimitedFiles

156 μs

using LinearAlgebra

135 μs

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.

md"""

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.

"""

129 μs

blurry_img

blurry_img = download("https://github.com/dbindel/cs4220-s20/raw/master/hw/code/proj1/data/blurry.png") |> load

944 ms

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).

md"""

## 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).

"""

152 μs

Matrix{RGB{N0f8}} (alias for Array{RGB{Normed{UInt8, 8}}, 2})

typeof(blurry_img)

5.7 μs

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.

md"""

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.

"""

155 μs

img_rgb_pixels (generic function with 1 method)

# Convenience function to count pixels

img_rgb_pixels(V) = prod(size(V)[2:3])

266 μs

arrays_to_img (generic function with 1 method)

# Convert 3-array representation to ordinary colorview representation

function arrays_to_img(V)

rgb = min.(max.(V, 0.0), 1.0)

return colorview(RGB, V[1,:,:], V[2,:,:], V[3,:,:])

end

397 μs

Now we can compute the 3-array representation of the image and illustrate the helpers.

md"""

Now we can compute the 3-array representation of the image and illustrate the helpers.

"""

121 μs

blurry_rgb

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.560784

# Compute the 3-array representation of the blurry image and illustrate the above helpers

blurry_rgb = channelview(float.(blurry_img))

38.6 ms

arrays_to_img(blurry_rgb)

938 μs

npixels = 79680

md"""npixels = $(img_rgb_pixels(blurry_rgb))"""

15.1 ms

Loading...