Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Project 2: Where are My Glasses?

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 $\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:

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?

![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
1.1 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

733 μs
using Images
5.8 s
using FileIO
212 μs
using FFTW
179 μs
using Colors
199 μs
using Roots
137 ms
using Plots
10.6 s
using DelimitedFiles
221 μs
using LinearAlgebra
151 μ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.
"""
145 μs
blurry_img
blurry_img = download("https://github.com/dbindel/cs4220-s20/raw/master/hw/code/proj1/data/blurry.png") |> load
1.8 s

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).
"""
186 μs
Matrix{RGB{N0f8}} (alias for Array{RGB{Normed{UInt8, 8}}, 2})
typeof(blurry_img)
28.9 μ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.
"""
161 μs
img_rgb_pixels (generic function with 1 method)
# Convenience function to count pixels
img_rgb_pixels(V) = prod(size(V)[2:3])
283 μ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
414 μ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.
"""
117 μ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))
33.0 ms
680 μs

npixels = 79680

md"""npixels = $(img_rgb_pixels(blurry_rgb))"""
14.6 ms
Loading...