DoF techniques

This week I started to research about DoF (Depth of Field) techniques, besides the classical gaussian blur I found 2 interesting techniques: the hexagonal DoF with bokeh from Frostbite2 and the Diffusion DoF used in Metro2033.
For the first, I was able to implement the hexagonal blur, but I yet didn't figured out how to compose the final image based on the blured image and the per-pixel CoC (Circle of Confusion). Simply interpolating between the original image and the blurred based on the CoC didn't looks good, and also there is the bleeding of in-focus/out-of-focus pixels colors. Here are some images (notice that I'm faking hdr to make brighter spots more perceptible)

Diffusion DOF

On the other hand the Diffusion DoF, in a elegantly way, uses the heat diffusion equation to implement the DoF effect. This method treat the image as a medium with each pixel color being its heat. Considering the pixel CoC as the heat conductivity of the pixel, we can simulate the heat diffusion  "spreading" the color among the pixels. Following we have the heat equation:

\frac{\partial u}{\partial t} = \beta \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)

Using the Alternate Direction Implicit (ADI) method we alternate the directions (x,y) solving two one dimensional problem at each time step. This 2D ADI method is unconditionally stable, this means that we can use any time step we want.
Simplifying the notation we have:

u_j^k = u_j^{k+1} + \frac{\beta \Delta t}{\Delta x^2}\( u_{j+1}^{k+1} - 2u_j^{k+1} + u_{j-1}^{k+1}\)

u_j^k = \(1+2s\)u_j^{k+1} - s\(u_{j+1}^{k+1} + u_{j-1}^{k+1}\)

where: s = \frac{\beta \Delta t}{\Delta x^2} , and u_{j}^{k} is the color in position j on time k.
 As the new value of u in the time k+1 depends on its adjacents values, we need to solve a system where:

\left( \begin{array}{cccccccc} b & c & & & & & & 0 \\ a & b & c & & & & & \\ & & & & & & & \\ & & & (...) & & & & \\ & & & & & & & \\ & & & & a & b & c & \\ & & & & & a & b & c\\ 0 & & & & & & a & b\\ \end{array} \right) \times \left( \begin{array}{c} u_0^{k+1} \\ u_1^{k+1} \\ \\ (...) \\ \\ \\ u_{n-1}^{k+1} \\ u_n^{k+1} \\ \end{array} \right) = \left( \begin{array}{c} u_0^k \\ u_1^{k} \\ \\ (...) \\ \\ \\ u_{n-1}^{k} \\ u_n^{k+1} \\ \end{array} \right)

with:

a = -s

b = 1 + 2s

c = -s

For this kind of matrix (tridiagonal)  we have fast methods to to solve the system, in my CPU implementation, I used the TDMA (also known as the Thomas algorithm).
Therefore we first have to calculate the diffusion in x direction and them in the y direction. Notice that for each row/column, we need to solve one system.

Before going to implement it in CrystalSpace, I implemented  it first in CPU to have a reference and here are the result:

Original image

CoC image

First step(diffusion in x direction)

Second step(diffusion in y direction)

You can download the code here.