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)
{
tbx_assert( grid.size() == mask.size());
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
// to 'mask'
int mask_nb_elts = mask.size().product();
std::vector<int> idx_list;
idx_list.reserve( mask_nb_elts );
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}} $$
Donate
No comments