*CMU 16-726 Spring 2021 Assignment #2*
Name: Juyong Kim
(#) Introduction
In this assignment, we implement Poisson blending, a gradient-domain processing technique that seamlessly blends an object or texture from a source image into a target image.
The insight behind Poisson blending is that people often care much more about the gradient of an image than the overall intensity.
So the method tries to change the source image so that the gradient of the source image is maximally preserved, while overall intensity is matched to the target image.
For the rest of the report, we first describe our method and implementation briefly.
Then we show the result of the toy problem given by the assignment.
The result of the blending including failure cases are shown following.
(#) Method
(##) Poisson Blending
Poisson blending method from the instruction is as follows.
Given the pixel intensities of the source image $\mathbf{s}$ and of the target image $\mathbf{t}$, we want to solve for new intensity values $\mathbf{v}$ within the source region $S$, which can be formulates as a least squares problem:
$$
\mathbf{v} = \arg \min_\mathbf{v} \sum_{i \in S, j \in N_i \cup S} \{(v_i - v_j) - (s_i - s_j)\}^2 + \sum_{i \in S, j \in N_i \cup \neg S} \{(v_i - t_j) - (s_i - s_j)\}
$$
where each $i$ is a pixel in the source region $S$, and $N_i$ is a 4-neighbor of $i$.
Each summation guides the gradient values to match those of the source region.
In the first summation, the gradients of the source image and the new image inside the region $S$ are compared.
In the second summation, the gradients of the source image and the gradient made from the new and target image on the boundary of $S$ are compared.
There are several considerations in running the algorithm.
First, the offset, or position, of the source image to be blended should be decided.
Although the offset is manually set, but there can be an algorithm to find the best offset so that the gradient at the boundary are similar between the source and the target image.
Second, for the blending output to be naturally looking, the texture, or the pattern, of the two images at the boundary should be similar.
Conceptually the method tries to preserve the gradient of the source image and the overall intensity of the source image can be changed towards the target image.
If the texture of two images differ significantly, there can be still noticable change at the boundary.
**Solving the Optimization**
The objective function above can be formularized as a least-square problem $\arg \min_\mathbf{v} \|A\mathbf{v}-b\|_2$, where each summand and the last row become a row.
However, the size of matrix $A$ is huge, even for the small image given by the assignment, solving it with the standard solution ($x = (A^\top A)^{-1}A^\top b$) takes long time.
Fortunately, since the matrix is very sparse (only 2 element per row is non-zero), we can solve with sparse matrix least-square solver.
Since we are dealing with RGB color images, we have to solve this optimization for each channel.
Here, the set $S$ is same for all channels, so $A$ is common and only $b$ is needed to be computed differently.
**Out-of-range Problem**
After solving the linear equation, some values of the solution may go out of the range of image intensity ($[0.0, 1.0]$ in this case).
We tried two methods to deal with the out-of-range values: Clipping and re-ranging.
For clipping, we can clip the pixel values greater than $1.0$ to $1.0$, and similarity for the values lower than $0.0$.
For re-ranging, we can apply linear point-processing so that the scale of each channel becomes $[0.0, 1.0]$.
After trying several examples, we found that clipping is enough when the blending is done smoothly, and the re-ranging is benefical in color2gray (described below).
(##) Toy Problem
To get a sense how to solve a least-square problem $\arg \min_x \|Ax-b\|_2$, we start with a toy problem of recovering an image with differences of gradients and a boundary condition.
The task here is that given an image $\mathbf{s}(x, y)$, we find the image $\mathbf{v}(x, y)$ that minimizes the following objective:
$$ \mathbf{v} = \arg\min_\mathbf{v} \bigg[ \sum_{x\in[H],y\in[W]} \{ (v(x+1,y)-v(x,y)) - (s(x+1,y)-s(x,y)) \}^2 \hspace{80pt} $$
$$ \hspace{40pt} + \sum_{x\in[H],y\in[W]} \{ (v(x+1,y)-t(x,y)) - (s(x+1,y)-s(x,y)) \}^2 + \{v(0,0) - s(0,0)\}^2 \bigg]. $$
This optimization finds an image which the gradients in both horizontal and vertical direction matches and the first pixel intensity.
We know the solution for $\mathbf{v}$ is the original image $\mathbf{s}$ making the objective function global minimum, zero.
(#) Bells and Whistles
(##) 1. Mixed Gradients
Same as the poisson blending, we can give different gradients than the source gradient.
Specifically, try the gradients of source or target image which has a larger absolute value.
So the objective function becomes as follows:
$$ \mathbf{v} = \arg \min_\mathbf{v} \sum_{i \in S, j \in N_i \cup S} \{(v_i - v_j) - d_{ij}\}^2 + \sum_{i \in S, j \in N_i \cup \neg S} \{(v_i - t_j) - d_{ij}\}, $$
$$ \text{where}~~~ d_{ij} = \begin{cases} s_i - s_j & ~~\text{if}~~ |s_i - s_j| \ge |t_i - t_j| \\ t_i - t_j & ~~\text{otherwise} \end{cases} $$
(##) 2. Color2Gray
One application of gradient-domain processing is to convert a color image into a gray-scale image that preserve the contrast between colors.
This is an example of tone-mapping problem, similar to that of converting HDR images to RGB displays.
To solve this, we first convert the RGB image into HSV color space, and apply mixed gradient method to the S (saturation) and V (value) channels.
Similarity to the toy problem, we add a squared difference between the first pixel of the output and the average of the S and V values of the first pixel.
(#) Results
(##) Part 1.1 Toy Problem
Figure [toy_input] and Figure [toy_output] show the input and the output of the toy problem.
The input is a small grayscale image with size $(110, 119)$.
But the number of constraints are proportional to the number of pixel and the matrix $A$ of the least square problem becomes $25952 \times 13090$, which is intractable to solve with standard method.
So the constraints are represented as sparse matrix `scipy.sparse.coo_matrix()` and the sparse matrix least-square solver `scipy.sparse.linalg.lsqr()` is used to solve the problem.
![Figure [toy_input]: Input](images/toy_problem.png width=200px) ![Figure [toy_output]: Output of `ssl.lsqr`