In [26]:
import numpy as np
import matplotlib.pyplot as plt
import skimage.io as io
import skimage as sk
In [27]:
from IPython.core.display import HTML
HTML("""
<style>

div.cell { /* Tunes the space between cells */
margin-top:1em;
margin-bottom:1em;
}

div.text_cell_render h1 { /* Main titles bigger, centered */
font-size: 2.2em;
line-height:0.9em;
}

div.text_cell_render h2 { /*  Parts names nearer from text */
margin-bottom: -0.4em;
}


div.text_cell_render { /* Customize text cells */
font-family: 'Georgia';
font-size:1.2em;
line-height:1.4em;
padding-left:3em;
padding-right:3em;
}

.output_png {
    display: table-cell;
    text-align: center;
    vertical-align: middle;
}

</style>

<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
The raw code for this IPython notebook is by default hidden for easier reading.
To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.

""")

#Trebuchet MS
Out[27]:
The raw code for this IPython notebook is by default hidden for easier reading. To toggle on/off the raw code, click here.

Final Project

Project1: High Dynamic Range

Modern cameras are unable to capture the full dynamic range of commonly encountered real-world scenes. In some scenes, even the best possible photograph will be partially under or over-exposed. Researchers and photographers commonly get around this limitation by combining information from multiple exposures of the same scene.

In this project, I implemented an algorithm that combines multiple exposures into a single high dynamic range radiance map. And then convert this randiance map to an image suiteable for display through tone mapping.

1.1 Radiance map construction

In this part, we firstly reover the radiance map from multiple exposures. The algorithm is based on Debevec and Malik, 1997. Here's the idea: we assume that the pixel values $Z_{ij}=f(E_i\Delta t_j)$, where i is a spatial index over pixels and j indexes over exposure times $\Delta t$. This equation can be rewriten as: $g(Z_{ij})=\ln E_i+\ln\Delta t_j$, where $g=\ln f^{-1}$. Then this problem can be formulated as finding $g(Z)$ that minimizes the following objective function: $\sum_i\sum_j{w(Z_{ij})[g(Z_{ij})-\ln E_i-\ln\Delta t_j]}^2+\lambda \sum_z[w(z) g(z)"]^2$.

The second term is a smoothness term on the sum of squared values of the second derivative of g to ensure that the function g is smooth; in this discrete setting we use $g(z)" = g(z-1)-2g(z)+g(z+1)$. This optimization problem can be further formulated as a least square problem. We also introduce a weighting function $w(z)$ to emphasize the smoothness and fitting terms toward the middle of the curve. $w(z)$ is chosen as a triangle function centered at 127.5.

By solving this optimization problem, we are able to find the function $g(Z)$. Here we show an exampled $g(Z)$:

In [28]:
im = io.imread("g(z).png")
f = plt.figure(figsize=(15,15))
ax = f.add_subplot(1,1,1)
ax.imshow(im)
plt.axis("off")
plt.show()

After we got the function g(Z), we can recover the radiance map E from it:

$\ln E_i=\frac{\sum_j w(Z_{ij})(g(Z_{ij})-\ln\Delta t_j)}{\sum_j w(Z_{ij})}$

We show some example radiance maps here:

In [29]:
im1 = io.imread("./radiance/arch_rad.png")
im2 = io.imread("./radiance/bonsai_rad.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
plt.axis("off")
plt.show()
In [30]:
im1 = io.imread("./radiance/chapel_rad.png")
im2 = io.imread("./radiance/house_rad.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
plt.axis("off")
plt.show()

Bell & Whistle -- image stack alignment

Some of the test image stacks are not well aligned. Therefore, I implemented an algorithm to automatically align the image stacks(images with different exposures) before running the HDR pipeline. The algorithm is based on multi-scale alignment similar to project 1 to search the translations along the image pyramid.

In [31]:
im1 = io.imread("./radiance/garden_rad_no_align.png")
im2 = io.imread("./radiance/garden_rad.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("Without alignment")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("With alignment")
plt.axis("off")
plt.show()

1.2 Tone Mapping

In this part, in order to recover an image at high dynamic range from the radiance map, we implemented two algorithms.

  1. Global Tone mapping: $V_{out}=\frac{V_{in}}{V_{in}+c}$, where $V_{in}$ is the radiance of the original pixel and $V_{out}$ is the brightness of the filtered pixel. This function will map the $V_{in}$ in the domain $[0,\infty)$ to an output range of $[0,1)$. c is used to control the brighness of the image.

  2. Local Tone Mapping based on Durand 2002. It uses a edge-preserving filter -- bilateral filter to display high-dynamic-range images. The algorithm is as follow:

(a) Your input is linear RGB values of radiance.

(b) Compute the intensity (I) by averaging the color channels.

(c) Compute the chrominance channels: (R/I, G/I, B/I)

(d) Compute the log intensity: L = log2(I)

(e) Filter that with a bilateral filter: B = bf(L)

(f) Compute the detail layer: D = L - B

(g) Apply an offset and a scale to the base: B' = (B - o) * s, where o = max(B), s = dR / (max(B) - min(B)).

(h) Reconstruct the log intensity: O = 2^(B' + D)

(i) Put back the colors: R',G',B' = O * (R/I, G/I, B/I)

(j) Apply gamma compression.

Here we show some example images with global(left) and local(right) Tone mappings.

In [32]:
im1 = io.imread("./global/arch_g.png")
im2 = io.imread("./local/arch_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("arch")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("arch")
plt.axis("off")
plt.show()
In [33]:
im1 = io.imread("./global/bonsai_g.png")
im2 = io.imread("./local/bonsai_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("bonsai")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("bonsai")
plt.axis("off")
plt.show()
In [34]:
im1 = io.imread("./global/chapel_g.png")
im2 = io.imread("./local/chapel_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("chapel")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("chapel")
plt.axis("off")
plt.show()
In [35]:
im1 = io.imread("./global/garage_g.png")
im2 = io.imread("./local/garage_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("garage")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("garage")
plt.axis("off")
plt.show()
In [36]:
im1 = io.imread("./global/garden_g.png")
im2 = io.imread("./local/garden_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("garden")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("garden")
plt.axis("off")
plt.show()
In [37]:
im1 = io.imread("./global/house_g.png")
im2 = io.imread("./local/house_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("house")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("house")
plt.axis("off")
plt.show()
In [38]:
im1 = io.imread("./global/mug_g.png")
im2 = io.imread("./local/mug_l.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("mug")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("mug")
plt.axis("off")
plt.show()

bilateral filtering

A bilateral filter is a non-linear, edge-preserving, and noise-reducing smoothing filter for images. It replaces the intensity of each pixel with a weighted average of intensity values from nearby pixels. This weight can be based on a Gaussian distribution.

In [39]:
im1 = io.imread("bf_before_filter.png")
im2 = io.imread("bf_after_filter.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,2,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("before bf filter")
ax2 = f.add_subplot(1,2,2)
ax2.imshow(im2)
ax2.set_title("after bf filter")
plt.axis("off")
plt.show()

project2: Image quilting

In this project, we implemented the image quilting algorithm for texture synthesis and transfer, described in SIGGRAPH 2001 by Efros and Freeman.

2.1: texture synthesis

In this part, we implemented an algorithm to create a larger texture image from a small sample. The main idea is to sample patches and lay them down in overlapping patterns, such that the overlapping regions are similar. The overlapping regions may not match exactly, which will result in noticeable edges. To fix this, I computed a path along pixels with similar intensities through the overlapping region and use it to select which overlapping patch from which to draw each pixel.

We compared the results from three approaches:

  1. Random Sampled Texture: it randomly samples square patches from the source image and generate an output image with a specific size. It starts from the upper-left corner, and tiles samples until the image is full.

  2. Overlapping Patches: It firstly samples a random patch for the upper-left corner. Then sample new patches to overlap with the existing one. Instead of randomly sampling, it searches all patches on the sample to find one with small difference with the overlapping region. It randomly choose on a patch whose cost is less than a threshold.

  3. Seam finding: In addition to the second approach, it removes the edge artifactes from the overlapping regions. We will explain this approach in the following part.

Illustration of seam finding

For two given patches sampled from the texture source, we firstly compute the cost (pixel-wised squared difference) between the overlapping regions. Then using dynamic programming, we can find the min-cost contiguous path from the uo to down (or left to right) side of the patch. This path defines a binary mask that specifies which pixels to copy from the newly sampled patch.

In [40]:
im1 = io.imread("illustration.png")
f = plt.figure(figsize=(15,15))
ax1 = f.add_subplot(1,1,1)
ax1.imshow(im1)
plt.axis("off")
plt.show()

Texture synthesis examples

We show some examples of texture synthesis for the above three approaches.

In [41]:
im1 = io.imread("./image_synthesis/apples.png")
im2 = io.imread("./image_synthesis/apples_random.png")
im3 = io.imread("./image_synthesis/apples_overlap.png")
im4 = io.imread("./image_synthesis/apples_cut.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax1.set_title("source", fontsize=20)
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2)
ax2.set_title("random sampling", fontsize=20)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
ax3.set_title("overlapping patches", fontsize=20)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
ax4.set_title("Seam finding", fontsize=20)
plt.axis("off")
plt.show()
In [42]:
im1 = io.imread("./image_synthesis/bricks_small.jpg")
im2 = io.imread("./image_synthesis/bricks_random.png")
im3 = io.imread("./image_synthesis/bricks_overlap.png")
im4 = io.imread("./image_synthesis/bricks_cut.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
plt.axis("off")
plt.show()
In [43]:
im1 = io.imread("./image_synthesis/cans.png")
im2 = io.imread("./image_synthesis/cans_random.png")
im3 = io.imread("./image_synthesis/cans_overlap.png")
im4 = io.imread("./image_synthesis/cans_cut.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
plt.axis("off")
plt.show()
In [44]:
im1 = io.imread("./image_synthesis/rice.png")
im2 = io.imread("./image_synthesis/rice_random.png")
im3 = io.imread("./image_synthesis/rice_overlap.png")
im4 = io.imread("./image_synthesis/rice_cut.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
plt.axis("off")
plt.show()
In [45]:
im1 = io.imread("./image_synthesis/lake.jpg")
im2 = io.imread("./image_synthesis/lake_random.png")
im3 = io.imread("./image_synthesis/lake_overlap.png")
im4 = io.imread("./image_synthesis/lake_cut.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
plt.axis("off")
plt.show()
In [46]:
im1 = io.imread("./image_synthesis/stone.jpg")
im2 = io.imread("./image_synthesis/stone_random.png")
im3 = io.imread("./image_synthesis/stone_overlap.png")
im4 = io.imread("./image_synthesis/stone_cut.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
plt.axis("off")
plt.show()

2.2 Texture Transfer

In this part, we implemented an algorithm for texture transfer based on texture synthesis approach. We create a texture sample that is guided by a pair of sample/target correspondence images. The image to be synthesized should satisfy two constraints:

  1. Each patch on the output image should come from the texture source image. Neighboring patches should have small cost in the overlapping region.
  2. The intensity difference betweeen the output image and the target image should be small.

We define the cost of each patch on the output image as:

$cost = \alpha*overlapping\_error + (1-\alpha)*correspondence\_map\_ssd$

The goal is to find patches with small costs within an error tolerance. Then we use above method to find the seam between two neighboring patches.

Bell&Whistle: iterative algorithm

I implemented the iterative texture transfer method described in the paper. Because of the added constraint, sometimes one synthesis pass through the image is not enough to produce a visually pleasing result. We iterate over the synthesized image several times, reducing the block size with each iteration. The only change from the non-iterative version is that in satisfying the local texture constraint the blocks are matched not just with their neighbor blocks on the overlap regions, but also with whatever was synthesized at this block in the previous iteration. We did $N=3$ iterations to produce the following results and reduced the patch size by a fourth each time and set $\alpha_i=\frac{i-1}{N-1}+0.1$.

We show some examples here:

In [47]:
im1 = io.imread("./texture_transfer/fabric.png")
im2 = io.imread("./texture_transfer/girl.png")
im3 = io.imread("./texture_transfer/fabric_girl.png")
im4 = io.imread("./texture_transfer/fabric_girl_iter.png")
f = plt.figure(figsize=(20,20))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
ax1.set_title("texture source", fontsize=20)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2, cmap="gray")
ax2.set_title("target", fontsize=20)
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
ax3.set_title("texture transfer", fontsize=20)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
ax4.set_title("iterative tecture transfer", fontsize=20)
plt.axis("off")
plt.show()
In [48]:
im1 = io.imread("./texture_transfer/cloth.jpg")
im2 = io.imread("./texture_transfer/einstein.jpg")
im3 = io.imread("./texture_transfer/cloth_einstein.png")
im4 = io.imread("./texture_transfer/cloth_einstein_iter.png")
f = plt.figure(figsize=(25,25))
ax1 = f.add_subplot(1,4,1)
ax1.imshow(im1)
plt.axis("off")
ax2 = f.add_subplot(1,4,2)
ax2.imshow(im2, cmap="gray")
plt.axis("off")
ax3 = f.add_subplot(1,4,3)
ax3.imshow(im3)
plt.axis("off")
ax4 = f.add_subplot(1,4,4)
ax4.imshow(im4)
plt.axis("off")
plt.show()

Difficulties of the project and what I learn:

I've worked hard on this project to improve the synthesized image quality. I found that it's related to many factors. For example, the choice of the parameters, e.g. tol, patch size, overlap size. It's also related to the source/target images I use. Sometimes, I need rescale the images so that the features' size between source and target are matched.

I learn a lot from this project and really have fun. It's a classic method for texture synthesis. I learned that even without any fancy method, e.g. neural network, we can still do a great job on image manipulation.

In [ ]: