2D biharmonic stencil a.k.a bilaplacian operator
Compute biharmonic weights on regular grids - 08/2014
Draft / work in progress
By applying the harmonic stencil to itself with the convolution operator we obtain the biharmonic stencil.
harmonic and biharmonic stencil (local host)
or here
For triangle mesh not sure where to look, perhaps here:
Biharmonic distance (based on cotan weights)
Localized bi-Laplacian Solver on a Triangle Mesh and Its Applications
Todo: 1d biharmonic stencil (1 -4 6 -4 1) (https://en.wikipedia.org/wiki/Finite_difference_coefficient) which you can infer by applying the harmonic stencil (1 -2 1) to itself
template< typename T>
void diffuse_biharmonic(Grid2_ref grid,
Grid2_const_ref mask,
int nb_iter = 128*16)
{
tbx_assert( grid.size() == mask.size());
Vec2_i neighs_1[4] = {Vec2_i( 1, 0),
Vec2_i( 0, 1),
Vec2_i(-1, 0),
Vec2_i( 0,-1)};
Vec2_i neighs_2[4] = {Vec2_i( -1, -1),
Vec2_i( -1, 1),
Vec2_i( 1, 1),
Vec2_i( 1, -1)};
// First we store every grid elements that need to be diffused according
// to 'mask'
int mask_nb_elts = mask.size().product();
std::vector 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 = T(0);
T mean = T(0);
for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs_1[i];
if( idx_neigh.is_in() /*check if inside grid*/ ){
T weight = (T)8;
sum = sum + grid( idx_neigh ) * weight;
mean += weight;
}
}
for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs_2[i];
if( idx_neigh.is_in() /*check if inside grid*/ ){
T weight = (T)-2;
sum = sum + grid( idx_neigh ) * weight;
mean += weight;
}
}
for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs_1[i] * 2;
if( idx_neigh.is_in() /*check if inside grid*/ ){
T weight = (T)-1;
sum = sum + grid( idx_neigh ) * weight;
mean += weight;
}
}
sum = sum / (T)( mean );
grid( idx ) = sum;
}
}
}
No comments