Image morphing with Delaunay triangulation

Question

i'm developing a simple algorithm to morph two images using keypoints and delaunay triangulation. The idea should be simple:

• choose the source control point
• choose the destination control point
• get the delaunay triangulation for source and destination frames
• for each pixel in the source image
  • get the pixel barycentric coordinates related to the source triangle in which the pixel lie
  • get the pixel barycentric coordinates related to the destination triangle in which the pixel lie
  • using the relationship Px = w1*v0x + w2*v1x + w3*v2x (the same for y, and the destination pixel) assign OUT[PdestX,PdestY] = IN[Px,Py].

But it doesn't work X_X This is my matlab source:

function out = myMorph(im1, p_source, p_dest, tri_source, tri_dest)

[h  w] = size(im1);

%get single column vectors for source and destination image control points
Psource_x   = p_source(:,1);
Psource_y   = p_source(:,2);
Pdest_x     = p_dest(:,1);
Pdest_y     = p_dest(:,2);

%for each intermediate frame...

out = zeros(size(im1));

%get triangles. Each array is 3n x 2, where n is the number of triangles
triangles_source = [];
triangles_dest = [];
for i= 1 : size(tri_source,1)
triangle_s = getTriangle(Psource_x,Psource_y,tri_source,i);
triangle_d = getTriangle(Pdest_x,Pdest_y,tri_dest,i);

triangles_source = cat(1,triangles_source,triangle_s);
triangles_dest = cat(1,triangles_dest,triangle_d);
end



%iterate each pixel
for x=1:h
for y=1:w

    %get the source and destination triangle for pixel [x y]

    %source triangle
    for t = 1 : 3 : size(triangles_source, 1)-2

       [w1,w2,w3,inTriangle] = inTri(x,y, ...
                                    triangles_source(t,1),triangles_source(t,2), ...
                                    triangles_source(t+1,1),triangles_source(t+1,2), ...
                                    triangles_source(t+2,1),triangles_source(t+2,2));
       if(inTriangle == 1)
           break;   %point [x,y] must belong to one (and only) triangle
       end
    end

    %source triangle
    for k = 1 : 3 : size(triangles_dest, 1)-2
       [w1d,w2d,w3d,inTriangleD] = inTri(x,y, ...
                                    triangles_dest(k,1),triangles_dest(k,2), ...
                                    triangles_dest(k+1,1),triangles_dest(k+1,2), ...
                                    triangles_dest(k+2,1),triangles_dest(k+2,2));
       if(inTriangleD == 1)
           break;
       end
    end

    v_source = [w1*triangles_source(t,1) + ...
                w2*triangles_source(t+1,1) + ...
                w3*triangles_source(t+2,1), ...
                w1*triangles_source(t,2) + ...
                w2*triangles_source(t+1,2) + ...
                w3*triangles_source(t+2,2)];

    v_dest = [w1d*triangles_dest(k,1) + ...
                w2d*triangles_dest(k+1,1) + ...
                w3d*triangles_dest(k+2,1),...
                w1d*triangles_dest(k,2) + ...
                w2d*triangles_dest(k+1,2) + ...
                w3d*triangles_dest(k+2,2)];


    if(inTriangle ~= 1 && inTriangleD ~= 1)
        continue;
    end

    v_source    = round(v_source);
    v_dest      = round(v_dest);

    if(v_source(1)>0 && v_source(1) <= h && ...
       v_source(2)>0 && v_source(2) <= w && ...
       v_dest(1)>0 && v_dest(1) <= h && ...
       v_dest(2)>0 && v_dest(2) <= w)

   disp('pixel warped')
    out(v_dest(1),v_dest(2)) = im1(v_source(1),v_source(2));
    end
   % else
    %    out(x,y) = im1(x,y); 

end
end

These are utility function for getting the control points

%Get control points used to morph im into another image
%im                     -> source image
%im2                    -> destination image
%linesNum               -> number of lines
function [P] = getControlPoints(im, controlPtsNum)

close all

P   = zeros(controlPtsNum, 2);

%select lines from source image
figure;
imshow(im,[]);title('select control points')

for i=1 : controlPtsNum
  %get source control point
 [x,y] =  ginput(1);
  P(i,:) = [x,y];

  hold on
   plot(x,y,'o','Color','r');
  hold off
end



%Get control points used to morph im into another image and do delaunay
%triangulation using the control points
%im                     -> source image
%im2                    -> destination image
%controlPtsNum          -> number of control points
function [P,tri] = getControlPointsAndTriangulate(im, controlPtsNum)

P = getControlPoints(im, controlPtsNum);

[h w] = size(im);

%Add corners to control points
P = cat(1, P, [1 1]);
P = cat(1, P, [w 1]);
P = cat(1, P, [1 h]);
P = cat(1, P, [w h]);

tri = delaunay(P(:,1),P(:,2));

hold on
triplot(tri,P(:,1),P(:,2))
hold on

and this function (i found on the net), test if a point lies on a given triangle, and return the u,v,w values:

function [w1,w2,w3,r] = inTri(vx, vy, v0x, v0y, v1x, v1y, v2x, v2y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inTri checks whether input points (vx, vy) are in a triangle whose
% vertices are (v0x, v0y), (v1x, v1y) and (v2x, v2y) and returns the linear
% combination weight, i.e., vx = w1*v0x + w2*v1x + w3*v2x and
% vy = w1*v0y + w2*v1y + w3*v2y. If a point is in the triangle, the
% corresponding r will be 1 and otherwise 0.
%
% This function accepts multiple point inputs, e.g., for two points (1,2),
% (20,30), vx = (1, 20) and vy = (2, 30). In this case, w1, w2, w3 and r will
% be vectors. The function only accepts the vertices of one triangle.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v0x = repmat(v0x, size(vx,1), size(vx,2));
v0y = repmat(v0y, size(vx,1), size(vx,2));
v1x = repmat(v1x, size(vx,1), size(vx,2));
v1y = repmat(v1y, size(vx,1), size(vx,2));
v2x = repmat(v2x, size(vx,1), size(vx,2));
v2y = repmat(v2y, size(vx,1), size(vx,2));
w1 = ((vx-v2x).*(v1y-v2y) - (vy-v2y).*(v1x-v2x))./...
((v0x-v2x).*(v1y-v2y) - (v0y-v2y).*(v1x-v2x)+eps);
w2 = ((vx-v2x).*(v0y-v2y) - (vy-v2y).*(v0x-v2x))./...
((v1x-v2x).*(v0y-v2y) - (v1y-v2y).*(v0x-v2x)+eps);
w3 = 1 - w1 - w2;
r = (w1>=0) & (w2>=0) & (w3>=0) & (w1<=1) & (w2<=1) & (w3<=1);  

Any suggestion? Bye!

Answer 1

I'm not able to reproduce the errors in your code because I don't have your input dataset, however, according to your description, you may have the same issue as I had when I attempted morphing an image by triangulation yesterday:

The number of triangles in the source triangulation and the destination triangulation are different.

This may be caused by what you described in your steps:

1. Perform Delaunay Triangulation with the source control points, get a triangular mesh
2. Perform Delaunay Triangulation with the destination control points, get a triangular mesh

Delaunay Triangulation is so clever that it uses minimal numbers of triangles for triangulation. It does not know the control points in step 2 are "transformed from" those in step 1. So the triangular meshes from steps 1 and 2 may contain different number of triangles! Here is an example and how to fix the problem:

Let's say you have constructed 2 lists of control points, "source CP" and "destination CP". "source CP" are the red points in case A. "destination CP" are the red points in cases B and C (they are identical).

Case A is obtained by executing Delaunay Triangulation over "source CP".

Case B is obtained by executing Delaunay Triangulation over "destination CP".

See? Case B contains 1 fewer triangle than case A!! If this happens you can't morph using the triangle lists in the Delaunay Triangulation of cases A and B.

A workaround is to get Case C with the same adjacency list and the same number of triangles as Case A, then you can perform image morphing with a pairwise triangle-to-triangle affine transformation approach.

Case C is obtained by just moving one control point in Case A, but keeping the same adjacency list.

Sure, the overlapping triangles have now become a new problem. I think you can place constraints like magnitude of distortion to prevent the triangles from overlapping. Besides, the triangle intersection test code you posted takes into account overlapping triangles by returning the triangle ID of the 1st triangle in the list that intersects with the query point.

So the point is you only need to execute Delaunay Triangulation once per source-destination transformation pair.

Hope this helps !