Matrix computations &
Numerical linear algebra

David Gleich

Purdue University

Fall 2012

Course number CS-51500

Tuesday and Thursday, 10:30-11:45am

Location Lawson 1106

Extra credit

$\newcommand{\mat}[1]{\boldsymbol{#1}} \renewcommand{\vec}[1]{\boldsymbol{\mathrm{#1}}} \newcommand{\vecalt}[1]{\boldsymbol{#1}} \newcommand{\conj}[1]{\overline{#1}} \newcommand{\normof}[1]{\|#1\|} \newcommand{\onormof}[2]{\|#1\|_{#2}} \newcommand{\itr}[2]{#1^{(#2)}} \newcommand{\itn}[1]{^{(#1)}} \newcommand{\eps}{\varepsilon} \newcommand{\kron}{\otimes} \DeclareMathOperator{\diag}{diag} \DeclareMathOperator{\trace}{trace} \DeclareMathOperator{\tvec}{vec} \newcommand{\prob}{\mathbb{P}} \newcommand{\probof}[1]{\prob\left\{ #1 \right\}} \newcommand{\pmat}[1]{\begin{pmatrix} #1 \end{pmatrix}} \newcommand{\bmat}[1]{\begin{bmatrix} #1 \end{bmatrix}} \newcommand{\spmat}[1]{\left(\begin{smallmatrix} #1 \end{smallmatrix}\right)} \newcommand{\sbmat}[1]{\left[\begin{smallmatrix} #1 \end{smallmatrix}\right]} \newcommand{\RR}{\mathbb{R}} \newcommand{\CC}{\mathbb{C}} \newcommand{\eye}{\mat{I}} \newcommand{\mA}{\mat{A}} \newcommand{\mB}{\mat{B}} \newcommand{\mC}{\mat{C}} \newcommand{\mD}{\mat{D}} \newcommand{\mE}{\mat{E}} \newcommand{\mF}{\mat{F}} \newcommand{\mG}{\mat{G}} \newcommand{\mH}{\mat{H}} \newcommand{\mI}{\mat{I}} \newcommand{\mJ}{\mat{J}} \newcommand{\mK}{\mat{K}} \newcommand{\mL}{\mat{L}} \newcommand{\mM}{\mat{M}} \newcommand{\mN}{\mat{N}} \newcommand{\mO}{\mat{O}} \newcommand{\mP}{\mat{P}} \newcommand{\mQ}{\mat{Q}} \newcommand{\mR}{\mat{R}} \newcommand{\mS}{\mat{S}} \newcommand{\mT}{\mat{T}} \newcommand{\mU}{\mat{U}} \newcommand{\mV}{\mat{V}} \newcommand{\mW}{\mat{W}} \newcommand{\mX}{\mat{X}} \newcommand{\mY}{\mat{Y}} \newcommand{\mZ}{\mat{Z}} \newcommand{\mLambda}{\mat{\Lambda}} \newcommand{\mPbar}{\bar{\mP}} \newcommand{\ones}{\vec{e}} \newcommand{\va}{\vec{a}} \newcommand{\vb}{\vec{b}} \newcommand{\vc}{\vec{c}} \newcommand{\vd}{\vec{d}} \newcommand{\ve}{\vec{e}} \newcommand{\vf}{\vec{f}} \newcommand{\vg}{\vec{g}} \newcommand{\vh}{\vec{h}} \newcommand{\vi}{\vec{i}} \newcommand{\vj}{\vec{j}} \newcommand{\vk}{\vec{k}} \newcommand{\vl}{\vec{l}} \newcommand{\vm}{\vec{l}} \newcommand{\vn}{\vec{n}} \newcommand{\vo}{\vec{o}} \newcommand{\vp}{\vec{p}} \newcommand{\vq}{\vec{q}} \newcommand{\vr}{\vec{r}} \newcommand{\vs}{\vec{s}} \newcommand{\vt}{\vec{t}} \newcommand{\vu}{\vec{u}} \newcommand{\vv}{\vec{v}} \newcommand{\vw}{\vec{w}} \newcommand{\vx}{\vec{x}} \newcommand{\vy}{\vec{y}} \newcommand{\vz}{\vec{z}} \newcommand{\vpi}{\vecalt{\pi}}$

Coloring within the lines with a matrix problem

Please work on the following problem and describe your work in complete sentences in a clearly prepared manuscript, submitted to me by December 7th. You must complete this work on your own.

In this assignment, we’ll implement the procedure described in Anat Levin, Dani Lischinski, and Yair Weiss in their paper “Colorization using optimization.” This paper describes how you can turn

into

by only specifying a few colors as in this picture

and solving two linear systems.

In this case, we’ll first see how some of the linear systems we’ve discussed in class arise in a real application. The challenge as an acting matrix analyst is, usually, to determine how to solve the problems faster. This is the task in the second question.

This assignment is heavily based on the authors’ own code package: http://www.cs.huji.ac.il/~yweiss/Colorization/ and my own work with Chris Maes on this problem. It’s a fun problem that’s perfect for an extra credit assignment as it builds off of the previous assignments in the class.

Problem 1: Quick demo! (5 points)

Solve a colorization problem yourself.

1. Download apple.mat, assemble.m, color_linsys_example.m You may need the function ntsc2rgb described below.
2. Run color_linsys_example (note, you need about 2GB of RAM for this).
3. Look at the source picture with imshow(src)
4. Can you tell if $\mA$ is ill-conditioned?

Problem 2: Generate the linear system (25 points)

In this first part of the extra-credit, we’ll build the linear system to solve. That is, you’ll write a function to take data from a file, and produce $\mA$, $\vb_1$ and $\vb_2$ for the above example.

Write the following function in Matlab

function [A,b1,b2] = color_linsys(src,marks,spread)
% COLOR_LINSYS
%
% Example:
%   src = double(imread('child.png')/255);
%   marks = double(imread('child-marked.png')/255);
%   [A,b1,b2] = get_color_linsys(src,marks);
%   x1 = A\b1;
%   x2 = A\b2;
%   final = assemble(src,x1,x2);
%   imshow(final);

We’ll walk through some details. One key to filling in the color is to separate the color information from the luminance information. This is how your TV used to transmit information. Given an RGB pixel value $\vp$ the YIQ color-space values are defined by the transformation:

$$\bmat{ y \\ i \\ q } = \bmat{ 1.0 0.956 0.621 \\ 1.0 -0.272 -0.647 \\ 1.0 -1.106 1.703 }^{-1} \bmat{ r \\ g \\ b }$$

After this transformation, the color information is contained in the values $i$ and $q$ – the information in $y$ only reflects brightness. This is why we only have to solve two linear systems, the value $\vb_1$ corresponds to information in the $i$ color channel and the value $\vb_2$ correspond to information in the $q$ color channel. Given the output from imread, we can convert it into yiq using the function rgb2ntsc available in the image processing toolbox. The core functionality of this function is really simple. The output from imread is an $m \times n \times 3$ array. (As an aside, I’m calling this an array rather than a 3-way matrix as it really just represents a bunch of data in memory.) Here are functions that will work instead:

function A = rgb2ntsc(A)
T = [1.0 0.956 0.621; 1.0 -0.272 -0.647; 1.0 -1.106 1.703];
[m, n] = size(A);
A = reshape(reshape(A, m*n, 3)\(T'), m, n, 3);

function A = ntsc2rgb(A)
T = [1.0 0.956 0.621; 1.0 -0.272 -0.647; 1.0 -1.106 1.703];
[m, n] = size(A);
A = reshape(reshape(A, m*n, 3)*T', m, n, 3);

So the first step of the function convert the image into ntsc (the same as yiq).

Now, we’ll build a pixel similarity map for the luminance data (the y channel). You can get this data by running

S = rgb2ntsc(src);
L = S(:,:,1);

Note that $\mL$ has the same dimensions as the image.
We’ll create a linear map that uses the luminance data to set the value of the pixel based on the surrounding pixels.

$$L_{i,j} = \sum_{i' = i-w}^{i+w} \sum_{j' = j-w}^{j+w} \exp(-(L_{i', j'}-L_{i,j})^2/(2\gamma\sigma_L^2) \text{excluding the $i,j$ term from the sum.}$$

The value $w$ is a “width” parameter that determines how far we’ll look for similarities in lumiance. We’ll set $w=1$ in this assignment so we only look at a $3\times 3$ patch around each pixel. The value $\sigma_L^2$ is the variance of luminance values in this $3\times 3$ patch of pixels including the pixel $L_{i,j}$.
Use the matlab function var to compute this score. For the picture above, I set $\gamma = 0.1$. The value of $\gamma$ determines how much color will spread throughout the picture. Making $\gamma$ smaller than $0.1$ is generally a bad idea because it’ll make your linear system ill-conditioned.

We’ve seen systems like this before in homeworks 1 and 2, so creating a matrix $\mA$ that implements this averaging function shouldn’t be hard. However, there is one more twist to this matrix. The above equation will diffuse color throughout the picture, but lose a little bit. We want to preserve color.
In order to do so, we need to ensure that each row of $\mA$ will sum to $0$ and that the diagonal entries of $\mA$ are all one. You can do this by rescaling the off-diagonals of the linear system. This is almost the matrix $\mA$ that you’ll return from the function. The difference is in how we handle the right hand side and the color information.

In order to do so, we need to determine which pixels have changed between the source image, which is black and white, and the marked image, which is in color. Here’s the function I used:

marked = sum(abs(S-M),3) > 1/10;

where $\mS$ and $\mM$ are the arrays for the source and marked image in the NTSC/YIQ color space, respectively. For each marked pixel, we want to change the linear system and treat it like a boundary term or source term. Again, we’ve seen this before! For each marked pixel, change the associated row in the matrix $\mA$ to assign that value to the right-hand side, i.e.

$$L_{i,j} = b_{i,j}$$

if the $i,j$ pixel is marked. For the right-hand side, we just use the color information from marked pixels. Thus, like explained above, you’ll get two right hand sides, one from M(:,:,2) and one from M(:,:,3).

To understand why this works, consider what the matrix $\mA$ does.
If you are a marked pixel, then it preserves your color information. If you aren’t a marked pixel, then your color is a weighted average of your neighbors color, where the weight is determine by the similarity of the luminance values. Since pixels are similar in regions that are all white, this will cause color to spread throughout the picture. If there are two marked regions, then they will spread until some mid-point.

Implement this function and make sure it correctly colors the child picture using the example code. Feel free to use child-res.png for reference.

Problem 3: Solving it fast (30 points)

Write the fastest solver you can to solve $\mA \vx = \vb$ with the two right-hand sides described above. You have a good deal of flexibility in your solution to this problem. For instance you could try various preconditioners for linear solvers, you could try different iterative methods, you could try finding the optimal omega for SOR, in short, treat this like a mini-report. That means it should be $3-5$ pages and use the tools you’ve acquired in this class to analyze the systems.

Give your runtime on the problems

• apple.png and apple-marked.png
• child.png and child-marked.png
• stanford.png and stanford-marked.png

If you wish to complete this portion independently of Problems 1 and 2, you may do so and contact me for the linear systems above. However, if I give them to you, then you may not complete Problems 1 or 2 for extra credit.

References

• Colorization using optimization
• Colorization using optimization (paper)
• Colorization using optimization (code)