CSE 559A: Computer Vision

Use Left/Right PgUp/PgDown to navigate slides

Fall 2018: T-R: 11:30-1pm @ Lopata 101

Instructor: Ayan Chakrabarti (ayan@wustl.edu).
Course Staff: Zhihao Xia, Charlie Wu, Han Liu

Sep 18, 2018

• Recitation will be this Friday (9/21) in Jolley 309.
• Will go over topics relevant to Pset.
• Problem Set due next Tuesday

# General

Notes on Parallelization

Say your simple code looked like this

for i in Indices:
y = func1(x[i])
z = func2(x[ifunc(i)])
p = func3(y,z)
out[i] = func4(p)

Replace as

y = func1(x[Indices])
z = func2(y[ifunc(Indices)])
p = func3(y,z)
out[Indices] = func4(p)
• Most numpy functions that act on single numbers can be used elementwise on arrays.
• Think about all the steps you would do for each number in a loop: Can these steps can be carried out independently for different loop indices ?
• If so, replace them with array operations.

# Image Restoration

$\hat{X} = \arg \min_X \|X - Y\|^2 + (X-\mu)^T \Sigma^{-1} (X-\mu)$

where, $$\mu = W^T\mu_c$$, $$\Sigma = W^TD_cW$$.

The solution is:

$\hat{X} = W^T(I+D_c^{-1})^{-1}W~~~(Y+W^TD_c^{-1}\mu_c)$

How would you code this up ?

# Image Restoration

$\hat{X} = W^T(I+D_c^{-1})^{-1}W~~~(Y+W^TD_c^{-1}\mu_c)$ $= W^T~~~{(I+D_c^{-1})^{-1} (WY + D_c^{-1}\mu_c)}$

# Image Restoration

$\hat{X} = W^T(I+D_c^{-1})^{-1}W~~~(Y+W^TD_c^{-1}\mu_c)$ $= W^T~~~{(I+D_c^{-1})^{-1} (\color{red}{WY} + D_c^{-1}\mu_c)}$

xc = wvlt(Y)



# Image Restoration

$\hat{X} = W^T(I+D_c^{-1})^{-1}W~~~(Y+W^TD_c^{-1}\mu_c)$ $= W^T~~~{(I+D_c^{-1})^{-1} \color{red}{(WY + D_c^{-1}\mu_c)}}$

xc = wvlt(Y)
xc += mu_c / sigma2c



# Image Restoration

$\hat{X} = W^T(I+D_c^{-1})^{-1}W~~~(Y+W^TD_c^{-1}\mu_c)$ $= W^T~~~\color{red}{(I+D_c^{-1})^{-1}}(WY + D_c^{-1}\mu_c)$

xc = wvlt(Y)
xc += mu_c / sigma2c
xc /= (1 + 1/sigma2c)



# Image Restoration

$\hat{X} = W^T(I+D_c^{-1})^{-1}W~~~(Y+W^TD_c^{-1}\mu_c)$ $= \color{red}{W^T}~~~{(I+D_c^{-1})^{-1}}(WY + D_c^{-1}\mu_c)$

xc = wvlt(Y)
xc += mu_c / sigma2c
xc /= (1 + 1/sigma2c)
X = invwvlt(xc)

Note that here mu_c and sigma2c are arrays the same size as xc.

# Image Restoration

More general optimization setting: $\hat{X} = \arg \min_X D(X; Y) + R(X)$

$\hat{X} = \arg \min_X \sum_n \|X[n]-Y[n]\|^2 + \lambda \sum_n \left(((G_x*X)[n])^2 + ((G_y*X)[n])^2\right)$

Let $$A_x$$ and $$A_y$$ be matrices corresponding to convolution with $$G_x$$ and $$G_y$$.

$\hat{X} = \arg \min_X \|X-Y\|^2 + \lambda~\left(\|A_xX\|^2 + \|A_yX\|^2\right)$

$=\arg \min_X (X-Y)^T(X-Y) + \lambda X^T(A_x^TA_x + A_y^TA_y)X$

Let's change this to our standard quadratic form:

$=\arg \min_X~~~X^T~\left(I + \lambda (A_x^TA_x + A_y^TA_y)\right)~X - 2X^TY + Y^TY$

$\mbox{And so,}~~~~~\hat{X} = \left(I + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~Y$

How do you do this matrix inverse ?

# Image Restoration

$\hat{X} = \left(I + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~Y$

Remember, (circular) convolution is diagonalized in the Fourier domain !

$A_x = S~D_x~S^*$

• Multiplication by $$S$$ is the inverse Fourier transform
• Multiplication by $$S^*$$ is the forward Fourier transform
• $$SS^* = I$$
• $$D_x$$ is a diagonal matrix with the Fourier transform coefficients of $$G_x$$

$A_x^TA_x = A_x^*A_x = S~D_x^*D_xS^* = S~\|D_x\|^2S^*$

$A_y^TA_y = S~\|D_y\|^2S^*$

$I+\lambda(A_x^TA_x+A_y^TA_y) = I + \lambda S (\|D_x\|^2+\|D_y\|^2) S^*$

$=S~\left(I + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)~S^*$

# Image Restoration

$\hat{X} = \left(I + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~Y$

$I+\lambda(A_x^TA_x+A_y^TA_y) = S~\left(I + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)~S^*$

$\left(I+\lambda(A_x^TA_x+A_y^TA_y)\right)^{-1} = S~\left(I + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)^{-1}~S^*$

$\hat{X} = S~\left(I + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)^{-1}~S^*Y$

Yf = fft2(Y)
Gxf, Gyf = fft2(Gx), fft2(Gy)
ratio = 1 + lambda*(np.abs(Gxf)**2 + np.abs(Gyf)**2)
Xf = Yf / ratio
X = ifft2(Xf)

What is this doing ?

It's down-weighting frequency components by a (real) factor where $$\|D_x\|^2$$ and $$\|D_y\|^2$$ are high.

Those are high for higher frequencies, because $$G_x$$ and $$G_y$$ are derivative filters.

So this operation down-weights higher frequency components $$\Rightarrow$$ Smooths the image.

# Image Restoration

Denoising

$\hat{X} = \arg \min_X \sum_n \|X[n]-Y[n]\|^2 + \lambda \sum_n \left(((G_x*X)[n])^2 + ((G_y*X)[n])^2\right)$

De-blurring

Say we know that our image has been blurred by a known blur kernel $$k$$

$Y[n] = (X*k)[n] + \epsilon[n]$

$\hat{X} = \arg \min_X \sum_n \|(X*k)[n]-Y[n]\|^2 + \lambda \sum_n \left(((G_x*X)[n])^2 + ((G_y*X)[n])^2\right)$

$\hat{X} = \arg \min_X \|A_kX-Y\|^2 + \lambda~\left(\|A_xX\|^2 + \|A_yX\|^2\right)$

where $$A_k$$ represents the action of convolution by blur kernel $$k$$.

$\hat{X} = \left(A_k^TA_k + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~A_k^TY$

Note that there is now $$A_k^TY$$ instead of just $$Y$$.

# Image Restoration

De-blurring / De-convolution

$\hat{X} = \arg \min_X \|A_kX-Y\|^2 + \lambda~\left(\|A_xX\|^2 + \|A_yX\|^2\right)$ $= \left(A_k^TA_k + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~A_k^TY$

We can do this in the Fourier domain again (assuming $$A_k$$ represents circular convolution).

$\hat{X} = S~\left(\|D_k\|^2 + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)^{-1}~S^*A_k^TY$

$\hat{X} = S~\left(\|D_k\|^2 + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)^{-1}~S^*~(SD_k^*S^*)Y$

$\hat{X} = S~\left(\|D_k\|^2 + \lambda \left(\|D_x\|^2+\|D_y\|^2\right)\right)^{-1}~D_k^*~S^*Y$

$X_f[u,v] = \frac{K_f[u,v]^*}{\|K_f[u,v]\|^2 + \lambda \left(\|G_{xf}[u,v]\|^2+\|G_{yf}[u,v]\|^2\right)} Y_f[u,v]$

If $$\lambda = 0$$, this reduces to:

$X_f[u,v] = \frac{K_f[u,v]^*}{\|K_f[u,v]\|^2} Y_f[u,v] = \frac{Y_f[u,v]}{K_f[u,v]},~\mbox{as}~\|K_f\|^2 = K_fK_f^*$

# Image Restoration

De-blurring / De-convolution

$\hat{X} = \arg \min_X \|A_kX-Y\|^2 + \lambda~\left(\|A_xX\|^2 + \|A_yX\|^2\right)$ $= \left(A_k^TA_k + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~A_k^TY$

$X_f[u,v] = \frac{K_f[u,v]^*}{\|K_f[u,v]\|^2 + \lambda \left(\|G_{xf}[u,v]\|^2+\|G_{yf}[u,v]\|^2\right)} Y_f[u,v]$

If $$\lambda = 0$$, this reduces to:

$X_f[u,v] = \frac{K_f[u,v]^*}{\|K_f[u,v]\|^2} Y_f[u,v] = \frac{Y_f[u,v]}{K_f[u,v]},~\mbox{as}~\|K_f\|^2 = K_fK_f^*$

But $$K_f[u,v]$$ may be zero or close to zero, in which case dividing will amplify noise.

So the Fourier transform of the kernel $$k$$ is telling us which frequency components
are severely attenuated by the kernel.

The "regularization" with $$\lambda > 0$$ helps stabilize the inversion, because if $$K_f[u,v]$$ is low for some $$(u,v)$$,
then the factor will downscale the input coefficient $$Y[u,v]$$.

This is called Wiener filtering or Wiener Deconvolution.

# Image Restoration

De-blurring / De-convolution

$\hat{X} = \arg \min_X \|A_kX-Y\|^2 + \lambda~\left(\|A_xX\|^2 + \|A_yX\|^2\right)$ $= \left(A_k^TA_k + \lambda (A_x^TA_x + A_y^TA_y)\right)^{-1}~~A_k^TY$

$X_f[u,v] = \frac{K_f[u,v]^*}{\|K_f[u,v]\|^2 + \lambda \left(\|G_{xf}[u,v]\|^2+\|G_{yf}[u,v]\|^2\right)} Y_f[u,v]$

If you're going to use this method to do de-convolution, you will have to account for the fact that your observed image was a "valid" convolution, and not a 'circular' convolution.

A good approximation is to do what's called "edge tapering". First pad the image $$Y$$ to a bigger image, where you make values go down smoothly to 0.

Then do the deconvolution, and crop out the central part.

# Image Restoration

We discussed cases when you know of a basis (wavelet / Fourier) where you can diagonalize your quadratic system matrix, and have a closed form expression for its inverse.

Not always the case. What if you wanted to exactly model valid convolution (not approximate it as circular) ? What if you observed values at a subset of pixels ?

Generally, what if you wanted to compute $$X = Q^{-1} Y$$ for some arbitrary symmetric positive-definite $$Q$$.

# Image Restoration

$X = Q^{-1}Y$

Let's consider a case where you can form $$Q$$.

• Never compute $$Q^{-1}$$, and then multiply by $$Y$$.
• Numerically unstable, more expensive.
• Call scipy.linalg.solve:
• Cholesky / LDL Decomposition: $$Q = L~D~L^T$$
• Always exists for a positive definite matrix. $$L$$ is lower triangular.
• Solve $$Qx = b \rightarrow LDL^Tx=b \rightarrow Ly = b, L^Tx = D^{-1}y$$

$\left[\begin{array}{ccccc}a & 0 & 0 & 0 & \ldots\\q&c&0&0&\ldots\\d&e&f&0&\ldots\\&&\vdots&&\end{array}\right]\left[\begin{array}{c}y_1\\y_2\\y_3\\\vdots\end{array}\right] = \left[\begin{array}{c}b_1\\b_2\\b_3\\\vdots\end{array}\right]$

# Image Restoration

$X = Q^{-1}Y$

More generally, when $$Q = A_1^TA_1 + A_2^TA_2 + \ldots$$, where $$A_1,A_2,\ldots$$ are sparse operations, that involve convolutions and element-wise operations.

• If $$A_1$$ is convolution with $$k$$, then you can get the effect of multiplying $$A_1^TA_1$$ with an image $$Y$$ by
• Convolving with $$k$$ first
• Then, convolving the result with the flipped version of $$k$$
• If $$A_1$$ is valid convolution, $$A_1^T$$ will correspond to "full" convolution with flipped version of $$k$$.
• If $$A_2$$ is convolution with $$k$$ followed by elmement-wise multiplication with a mask image, then $$A_2^TA_2$$ is
• Convolution with $$k$$
• Convolution with flipped version of $$k$$
• So even when we can't form $$Q$$, we can carry out the actions $$QY$$, as well as $$Z^TQY$$
• Compute $$QY$$
• Take element-wise product of the result with $$Z$$ and sum.

# Image Restoration

$x = Q^{-1}b$

Solve by the Conjugate Gradient method.

• Generic algorithm for solving $$Qx = b$$ for symmetric positive definite $$Q$$.
• Useful when you can multiply by $$Q$$ but not 'form' it.

Basic Idea

• For a given set of vectors $$\{p_1,p_2,\ldots p_N\}$$
• that are same size as $$x$$
• linearly independent
• $$N$$ = dimensionality of $$x$$
• We can write any $$x = \sum_i \alpha_i p_i$$
• If we also choose the vectors to be 'conjugate' such that $$p_i^TQp_j = 0$$ for $$i\neq j$$:

$Qx = b \rightarrow p_k^TQx = p_k^Tb \rightarrow \alpha_i p_k^TQp_k = p_k^Tb \rightarrow \alpha_i = \frac{p_k^Tb}{p_k^TQp_k}$

Iterative Algorithm

• Begin with some guess $$x_0$$ for $$x$$ (say all zeros)
• $$k = 0, r_0 \leftarrow b-Qx_0,~~p_0\leftarrow r_0$$
• Repeat
• $$\alpha_k \leftarrow \frac{r_k^Tr_k}{p_k^TQp_k}$$
• $$x_{k+1} = x_k + \alpha_k p_k$$
• $$r_{k+1} = r_k - \alpha_k Qp_k$$
• $$\beta_k = \frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}$$
• $$p_{k+1} = r_{k+1} + \beta_kp_k$$
• $$k=k+1$$

Think about what you would do when: $$Q =\left(A_k^TA_k + \lambda (A_x^TA_x + A_y^TA_y)\right),~~~~b=A_k^TY$$