Laplacian interpolation implementation in matlab

Some time ago, I quoted a passage from a website which is able to give a general idea about the Laplacian interpolation.

The address is http://blog.csdn.net/seamanj/article/details/50701112

Here, I want to implement this algorithm in matlab.

Firstly, I would like to elaborate the underlying mathematical principle.

The requirement originates from a paper named “Anatomy transfer“, the 7th section of which tells us: Given by value on the boundary, we can use Laplace interpolation to infer the value in the interior part.



Then I recommend to refer to the 2.1 part of alec jacobson thesis. In order to better understand, we reinterpret this part in our own way showed in N4 of Alec jacobson thesis analysis. From which we can get the following formula:

LΩ¯ΩuΩ=MΩ¯Ω¯fΩ¯−LΩ¯0u0

where Ω¯¯¯=Ω+∂Ω=Ω+0

Here, our purpose is to make Δu zero; therefore, fΩ¯=0, which leads to the following formula:

LΩ¯ΩuΩ=−LΩ¯0u0

where LΩ¯Ω,LΩ¯0 are parts of LΩ¯Ω¯, which is a stiffness matrix being able to be solved through solve stiffness matrix in matlab, u0 is the known value on the boundary. That leaves only uΩ as the unknown.

the main.m file is listed as follows, the complete source code can be found in attachment.

xRes = 16;
yRes = 16;

wrap = 0;
[F,V,res,edge_norms] = create_regular_grid(xRes,yRes,wrap,wrap);
V = [V(:,1), V(:,2), zeros(size(V,1),1)];
subplot_handle = subplot(2,2,1);
cla(subplot_handle);
tsurf(F,V,0,0);
title('Input domain');

indices = 1:size(V,1);
border = indices( V(:,1) == 0 | V(:,1) == 1 | V(:,2) == 0 | V(:,2) == 1 );
interior = indices(~ismember(indices, border));

subplot_handle = subplot(2,2,2);
cla(subplot_handle);
display_domain(F,V,border);
title('Paritioned domain');

[bi_L,bi_U,bi_P,bi_Q,bi_R,bi_S] = biharm_factor_system( ...
V, ...
F, ...
interior, ...
border);

V(border,:) = [ V(border,1).*2 ...
V(border,2).*2 ...
-2 * V(border,1) .*V(border,1) + 2 * V(border,1) ];
%zeros(size(border,2),1)

% Build RHS and use factor system matrix to solve for new positions in each
% coordinate
bi_V = biharm_solve_with_factor( ...
bi_L, bi_U, bi_P, bi_Q, bi_R, bi_S, ...
V, interior, border);
%
% % Display solution
subplot_handle = subplot(2,2,3);
% clear cubplot just in case
cla(subplot_handle);
display_domain(F,bi_V,border);
view(3)
zoom out
zoom(1.5)
title('laplacian interpolation');
axis equal;
axis vis3d;

The result picture is like:



we only change the positions of the border vertices expressed in dark red dot, the interior part in yellow color would be optimized by laplacian interpolation.

the complete matlab source code