# 2D harmonic stencil a.k.a Laplace operator

Related notes on biharmonic weights for grids

1D harmonic stencil:

(1 -2 1)


The 2D harmonic stencil (a.k.a the Laplace operator) is:

$\Delta f(x,y) \approx \frac{f(x-h,y) + f(x+h,y) + f(x,y-h) + f(x,y+h) - 4f(x,y)}{h^2},$

1

1  -4   1

1


Diffusion algorithm on a 2D grid that will converge to harmonic weights:

/// @param grid : grid of values/vectors you want to diffuse
/// @param mask : values to '0' are not diffused
/// @tparam T : must implement constructor T( double ) and overload
/// 'T operator/(double)' 'T operator+(T)'
/// @note: diffusing is equivalent to solving Lx=b using Gauss Seidel iterations.
/// (L is the laplacian of our grid)
template< typename T>
void diffuse_laplacian(Grid2_ref<T> grid, Grid2_const_ref<int> mask, int nb_iter = 128*16)
{
Vec2_i neighs[4] = {Vec2_i( 1, 0),
Vec2_i( 0, 1),
Vec2_i(-1, 0),
Vec2_i( 0,-1)};

// First we store every grid elements that need to be diffused according
std::vector<int> idx_list;
for(Idx2 idx(mask.size(), 0); idx.is_in(); ++idx) {
if( mask( idx ) == 1 )
idx_list.push_back( idx.to_linear() );
}

// Diffuse values with a Gauss-seidel like scheme
for(int j = 0; j < nb_iter; ++j)
{
// Look up grid elements to be diffused
for(unsigned i = 0; i < idx_list.size(); ++i)
{
Idx2 idx(grid.size(), idx_list[i]);

T sum(0.);
int nb_neigh = 0;

for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs[i];
if( idx_neigh.is_in() /*check if inside grid*/ ){
sum = sum + grid( idx_neigh );
nb_neigh++;
}
}

sum = sum / (T)( nb_neigh );
grid( idx ) = sum;
}
}
}


Recall that one iteration of Guauss seidel to solve a linear system $\mathbf {Ax} = \mathbf{b}$ is:

$$x_i^{k+1}=\frac{ \overbrace{\sum \limits_{j=1}^{i-1} {a_{ij}.x_j^{k+1}}}^{ \text{already computed value} \\ \text{in previous iteration} }+ (b_i-\sum \limits_{j=i+1}^{n}a_{ij}.x_j^k) }{a_{ii}}$$

(optional field, I won't disclose or spam but it's necessary to notify you if I respond to your comment)
All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Anti-spam question: