2D Quadrilateral Interpolation: Barycentric Coordinates

Contents

Interpolation on a general quadrilateral

When dealing with general quadrilaterals, the explicit expression of the shape functions is not easy to found. But we have the possibility of computing the value of such functions on a point p using the barycentric coordinates for quadrilaterals.

Similar to the triangular interpolation, now consider a point $$ p=(x,y)$, in a quadrilateral defined by four vertices $$v_{j}$ numbered counterclokwise as usual(see figure). One can define the barycentric coordinades of this point which are equivalent to the value of the 2D shape functions $$\psi_i^k(x,y)$ defined for this quadrilateral $$Q^k$:

satisfying $$\psi_i^k(v_j)= \delta_{ij}$.

This way, one can express any point using four coordinates (barycentric) $$p=(\alpha_1,\alpha_2,\alpha_3,\alpha_4)$. In particular $$v_1=(1,0,0,0); v_2=(0,1,0,0); v_3=(0,0,1,0); v_4=(0,0,0,1)$.

The barycentric coordinates satisfy $$\alpha_1+\alpha_2+\alpha_3+\alpha_4=1$. Moreover all $$\alpha_i \geq 0$ unless the point is outside of the quadrilater.

Thus, the coordinates of the point

$$ p=\alpha_1 v_1+\alpha_2 v_2+\alpha_3 v_3+\alpha_4 v_4$.

To compute the values of $$\alpha_j$ we re-write the previous equation as

$$ p=(1-\lambda)(1-\mu)v_1+\lambda(1-\mu) v_2+\lambda\mu v_3+(1-\lambda)\mu v_4$.

where $$\lambda, \mu \in [0,1]$

Therefore,

$$\alpha_1= (1-\lambda)(1-\mu), \quad \alpha_2= \lambda(1-\mu), \quad \alpha_3= \lambda\mu, \quad \alpha_4= (1-\lambda)\mu$

Now to compute these two parameters, we rename the terms in the previous equation as:

$$\mathbf{a}+\mathbf{b}\lambda+\mathbf{c}\mu+\mathbf{d}\lambda\mu=0$

where $$\mathbf{a}= v_1-p, \quad \mathbf{b}= v_2-v_1, \quad \mathbf{c}= v_4-v_1, \quad \mathbf{d}= v_1-v_2-v_3+v_4.$

Finally we make a cross product of the last equation with $$\mathbf{a}+\mathbf{c}\mu$, getting

$$(\mathbf{a}+\mathbf{c}\mu) \times (\mathbf{b}+\mathbf{d}\mu) = 0$

or

$$(\mathbf{c} \times \mathbf{d})\mu^2+(\mathbf{c} \times \mathbf{b}+\mathbf{a} \times \mathbf{d})\mu+\mathbf{a} \times \mathbf{b}=0$

which can be solved (Hint. The case $$\mathbf{c} \times \mathbf{d} =0$ must be treated separately)

The equation for $$\lambda$ is obtained analogously

$$(\mathbf{b} \times \mathbf{d})\lambda^2+(\mathbf{b} \times \mathbf{c}+\mathbf{a} \times \mathbf{d})\lambda+\mathbf{a} \times \mathbf{c}=0$

Compute Barycentric Coordinates on a general quadrilateral

Given a quadrilateral defined by vertices $$v_1=(0,0), v_2=(5,-1), v_3=(4,5), v_4=(1,4)$, compute the barycentric coordinates of point $$p=(3,2)$.

Verify that $$ p=\alpha_1 v_1+\alpha_2 v_2+\alpha_3 v_3+\alpha_4 v_4$.

% Data values
p=[3 2];
v1=[0 0];
v2=[5 -1];
v3=[4 5];
v4=[1 4];
% define the previous notation
a=(v1-p);
b=v2-v1;
c=v4-v1;
d=v1-v2-v4+v3;

% compute 2on order equation A mu^2 + B mu + C=0
% as the vertices are 2D, we add a zero third component
% to compute cross products.
A=cross([c,0],[d,0]); %must be 3D vectors
B=cross([c,0],[b,0])+cross([a,0],[d,0]);
C=cross([a,0],[b,0]);
% Only third component is needed (the other two are zero)
 A=A(3);
 B=B(3);
 C=C(3);
%
% Check for unique solutions
%
if (abs(A)<1.e-14)
    u1= -C/B;
    u2=u1;
else
    %
    % Check for non complex solutions
    %
    if (B^2-4*A*C >0)
        u1=(-B+sqrt(B^2-4*A*C))/(2*A);
        u2=(-B-sqrt(B^2-4*A*C))/(2*A);
    else %complex solution
        u1=-1000;
        u2=u1;
    end
end
%
mu=-10000; %stupid value small enough
if(u1>=0 && u1<=1)
    mu=u1;
end
if(u2>=0 && u2<=1)
    mu=u2;
end

% compute 2on order equation A lambda^2 + B lambda + C=0
A=cross([b,0],[d,0]); %must be 3D vectors
B=cross([b,0],[c,0])+cross([a,0],[d,0]);
C=cross([a,0],[c,0]);
% Only third component is needed (the other two are zero)
A=A(3);
B=B(3);
C=C(3);
%
% Check for unique solutions
%
if (abs(A)<1.e-14)
    w1= -C/B;
    w2=w1;
else
    %
    % Check for non complex solutions
    %
    if (B^2-4*A*C >0)
        w1=(-B+sqrt(B^2-4*A*C))/(2*A);
        w2=(-B-sqrt(B^2-4*A*C))/(2*A);
    else %complex solution
        w1=-1000;
        w2=w1;
    end
end
%
lambda=-10000; %stupid value
if(w1>=0 && w1<=1)
    lambda=w1;
end
if(w2>=0 && w2<=1)
    lambda=w2;
end
[mu,lambda] %parameters
% Barycentric coordinates
alpha1=(1-mu)*(1-lambda);
alpha2=lambda*(1-mu);
alpha3=mu*lambda;
alpha4=(1-lambda)*mu;
alphas=[alpha1,alpha2,alpha3,alpha4]
% obtained point
q=alpha1*v1+alpha2*v2+alpha3*v3+alpha4*v4;
p-q %we must recover the same point
% Plot the results
vertices=[v1,v2,v3,v4];
plotRectangle(v1,v2,v3,v4);
hold on;
plot(p(:,1),p(:,2),'o');
plot(q(:,1),q(:,2),'*');
hold off;
ans =

   5.0000e-01   6.2500e-01


alphas =

   1.8750e-01   3.1250e-01   3.1250e-01   1.8750e-01


ans =

     0     0

Exercise 1:

From the previous code implement the function baryCoordQuad (the Quadrilateral version of the baryCoord function).

function [alphas,isInside]=baryCoordQuad(vertices,point)

Exercise 2:

From the mesh file meshTwoHolesQuad.m, find the rectangle where the point p=[39,7] belongs to.

Sol:

i = 404

vertices =

41.9240    8.7481
37.5060    8.6827
36.5394    4.7543
41.0196    4.8224

(c)Numerical Factory 2018