Generalized Barycentric Coordinates for Warpping¶

$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This tours tests several barycentric coordinates (mean value, harmonic and green) for non-convex polygons, and apply them to 2D warping of images.

In [2]:

addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/meshdeform_4_barycentric')

Domain to Warp¶

Create a cage that bound the domain to warp.

In [3]:

cage = 'V';

First we create a 2D closes polygon, which will be a cage used to perform 2D shape deformation.

In [4]:

delta = .03;
rho = .2;
eta = .05;
x1 = .5-rho-eta/2; x2 = .5-eta/2;
x3 = .5+eta/2; x4 = .5+rho+eta/2;
switch cage
    case 'L'
        V = [[delta;delta] [1-delta;delta] [1-delta;delta+rho] [delta+rho;delta+rho] [delta+rho;1-delta] [delta;1-delta]];
    case 'U'
        V = [[x1;delta] [x4;delta] [x4;1-delta] [x3;1-delta] ...
            [x3;delta+rho] [x2;delta+rho] [x2;1-delta] [x1;1-delta] ];
    case 'V'
        V = [[x1;delta] [x4;delta] [x4;1-delta] [x3;1-delta] ...
            [.5;delta+rho] [x2;1-delta] [x1;1-delta] ];
end
k = size(V,2);

Compute a grid.

In [5]:

n = 200;
x = linspace(0,1,n);
[Y,X] = meshgrid(x,x);

Indicator of the shape.

In [6]:

S = 1 - inpolygon(X,Y,V(1,:),V(2,:));

Compute a check-board texture inside the L shaped domain.

In [7]:

m  = 30;
[XT,YT] = meshgrid((0:n-1)/n,(0:n-1)/n);
T = mod( floor(XT*m)+floor(YT*m),2 );
T(S==1) = 1;

Display it.

In [8]:

lw = 3; ms = 25;
clf; hold on;
plot_surf_texture(cat(3,X,Y,zeros(n)), T');
h = plot(V(1,[1:end 1]), V(2,[1:end 1]), 'r.-');
set(h, 'LineWidth', lw); set(h, 'MarkerSize', ms);
view(2); axis('off'); axis('equal');

Mean-valued Coordinates¶

Mean valued coordinates are usually applied for 3D mesh parameterization, and are used to approximate the Lapalcian on a 1-ring of a triangulation.

The mean value coordinates for star-shaped polygons where introduced in

M. S. Floater, Mean value coordinates Comp. Aided Geom. Design 20, 19-27, 2003.

There extension to arbitrary polygon is presented in

K. Hormann and M.S. Floater, Mean value coordinates for arbitrary planar polygons ACM Transactions on Graphics, 25 p. 1424-1441, 2006.

See also the following paper for application to shape warping:

T. Ju, S. Schaefer and J. Warren, Mean Value Coordinates for Closed Triangular Meshes ACM SIGGRAPH 2005, pages 561-566, 2005.

Useful operator.

In [9]:

dotp = @(a,b)sum(a.*b);
crossp = @(a,b)a(1,:).*b(2,:)-a(2,:).*b(1,:);
normalize = @(a)a./repmat(sqrt(sum(a.^2)), [2 1]);

Points of the domain.

In [10]:

W = [X(:)';Y(:)'];

Each C(:,:,i) will be a barycentric coordinate weight.

In [11]:

C = zeros(n,n,k);

Select a point on the polygon.

In [12]:

i = 4;
vi = V(:,i);

Compute the mean coordinate at each point location using the tangent of the two adjacent angles.

In [13]:

U = repmat(vi,[1 n^2])-W;
nb = normalize( U );

length

In [14]:

d = sqrt( sum(U.^2) );
for j=mod([i-2,i],k)+1
    % point
    vj = V(:,j);
    na = normalize( repmat(vj,[1 n^2])-W );
    % angle
    dp = dotp(na,nb);
    theta = acos(clamp(dp,-1,1));
    % add tangent of half angle
    C(:,:,i) = C(:,:,i) + reshape( tan(theta/2) ./ d, [n n]);
end

Exercise 1

Compute the full set of mean coordinates.

In [15]:

exo1()

In [16]:

%% Insert your code here.

Normalize them.

In [17]:

C = C ./ repmat( sum(C,3), [1 1 k] );

Extract one of the coordinates, set it to zero outside.

In [18]:

i = 4;
c = abs(C(:,:,i))+1e-3;
c(S==1) = 0;

Display it.

In [19]:

nl = 15;
t = linspace(0,1,n);
B = display_shape_function(c');
clf; hold on;
imagesc(t,t,B); axis('image'); axis('off');
contour(t,t,c',nl, 'k');
colormap jet(256);
h = plot(V(1,[1:end 1]), V(2,[1:end 1]), 'r.-');
set(h, 'LineWidth', lw); set(h, 'MarkerSize', ms);

Interpolation function.

In [20]:

applyinterp = @(C,x)sum(repmat(reshape(x(:), [1 1 k]),[n n 1]).*C,3);

Apply the interpolation weight to the X/Y coordinate to test for the linear precision of the coordinates.

In [21]:

clf;
imageplot(applyinterp(C,V(1,:)), 'Should be X', 1,2,1);
imageplot(applyinterp(C,V(2,:)), 'Should be Y', 1,2,2);
colormap jet(256);

Final position of the cage.

In [22]:

V2 = V;
switch cage,
    case 'L'
        V2(:,4) = [1-delta;1-delta];
    case {'U' 'V'}
        V2(:,3) = [1-delta;delta];
        V2(:,4) = [1-delta;delta+rho];
end

Modify the position of the cage.

In [23]:

rho = .7;
V1 = V*(1-rho) + V2*rho;

Warp the grid.

In [24]:

X1 = applyinterp(C,V1(1,:));
Y1 = applyinterp(C,V1(2,:));

Display the warped texture.

In [25]:

clf; hold on;
plot_surf_texture(cat(3,X1,Y1,zeros(n)), T');
h = plot(V1(1,[1:end 1]), V1(2,[1:end 1]), 'r.-');
set(h, 'LineWidth', lw); set(h, 'MarkerSize', ms);
view(2); axis('off'); axis('equal');
axis([0,1,0,1]);

Harmonic Coordinates¶

Mean valued coordinates goes in some sense "through" the cage. To prevent this, one can compute coordinates using a diffusion within the cage. This leads to harmonic coordinates.

The harmonic coordinates are introduced in

P. Joshi, M. Meyer, T. DeRose, B. Green and T. Sanocki Harmonic coordinates for character articulation ACM Trans. Graph, 3(26), p.71, 2007

Compute the points on the grid that are along the boundary, together with the curvilinear absice.

In [26]:

Vi = round(V*n);
bound = Vi(:,1);
loc = 1;
abscur = 0;
for i=1:k
    j = mod(i,k)+1;
    d = round(norm(Vi(:,i)-Vi(:,j)));
    t = repmat((1:d)/d, [2 1]);
    bound = [bound repmat(Vi(:,i),[1 d]).*(1-t) + repmat(Vi(:,j),[1 d]).*t];
    abscur = [abscur abscur(end)+(1:d)/d];
    loc(end+1) = loc(end)+d;
end
q = size(bound,2);
bound = round(bound);
I = bound(1,:) + (bound(2,:)-1)*n;

Initialize the coordinates.

In [27]:

C = zeros(n,n,k);

Compute the "hat" interpolation function at a given vertex location.

In [28]:

i = 4;
u = zeros(k+1,1);
u(i) = 1;
if i==1
    u(end)=1;
end
u = interp1(0:k,u, abscur);

Initialize the coordinate.

In [29]:

Ci = zeros(n);

Perform diffusion.

In [30]:

sel1 = [2:n 1];
sel2 = [n 1:n-1];
Ci = ( Ci(sel1,:) + Ci(:,sel1) + Ci(sel2,:) + Ci(:,sel2) )/4;

Impose value.

In [31]:

Ci(I) = u;

Display the first iteration.

In [32]:

clf;
imageplot(Ci');
colormap jet(256);

Exercise 2

Perform the full computation of the coordinate |C(:,:,i)| by iterating the diffusion and imposing the boundary value.

In [33]:

exo2()

In [34]:

%% Insert your code here.

Display one the functions.

In [35]:

c = C(:,:,i); c(S==1) = 0;
nl = 15;
t = linspace(0,1,n);
B = display_shape_function(c');
clf; hold on;
imagesc(t,t,B); axis('image'); axis('off');
contour(t,t,c',nl, 'k');
colormap jet(256);
h = plot(V(1,[1:end 1]), V(2,[1:end 1]), 'r.-');
set(h, 'LineWidth', lw); set(h, 'MarkerSize', ms);

Exercise 3

Compute the full set of coordinate functions |C|.

In [36]:

exo3()

In [37]:

%% Insert your code here.

Apply the interpolation weight to the X/Y coordinate to test for the linear precision of the coordinates. Here it is only valid inside the shape.

In [38]:

clf;
A = applyinterp(C,V(1,:)); A(S==1) = 0;
imageplot(A, 'Should be X', 1,2,1);
A = applyinterp(C,V(2,:)); A(S==1) = 0;
imageplot(A, 'Should be Y', 1,2,2);
colormap jet(256);

Modify the position of the cage.

In [39]:

rho = .7;
V1 = V*(1-rho) + V2*rho;

Warp the grid.

In [40]:

X1 = applyinterp(C,V1(1,:));
Y1 = applyinterp(C,V1(2,:));
X1(S==1) = Inf; Y1(S==1) = Inf;

Display the warped texture.

In [41]:

clf; hold on;
plot_surf_texture(cat(3,X1,Y1,zeros(n)), T');
h = plot(V1(1,[1:end 1]), V1(2,[1:end 1]), 'r.-');
set(h, 'LineWidth', lw); set(h, 'MarkerSize', ms);
view(2); axis('off'); axis('equal');
axis([0,1,0,1]);

Green Coordinates¶

Another set of coordinates, that are fast to compute and does not go "though" the cage, are the Green Coordinates, introduced in

Y. Lipman, D. Levin and D. Cohen-Or, Green Coordinates ACM Trans. Graph., 27(3), pages 1-10, 2008}.

Fist compute the oriented normal to the cage.

In [42]:

N = V(:,[2:end 1]) - V;
N = N ./ repmat( sqrt(sum(N.^2)), [2 1] );
N = -[-N(2,:); N(1,:)];

Display the normals.

In [43]:

clf; hold on;
h = plot(V(1,[1:end 1]), V(2,[1:end 1]), 'r.-');
set(h, 'LineWidth', lw); set(h, 'MarkerSize', ms);
rho = .1;
for i=1:k
    j = mod(i,k)+1;
    a = mean( V(:,[i,j]), 2);
    h = plot( [a(1) a(1)+rho*N(1,i)], 