Face Morphing with Mathematica

Question

Can Mathematica be used to produce animations of face morphing (as described in wikipedia?

This is a self-answered question because I felt there was a lack of resources on this topic, to share funny stuff and benefit from experienced users' advice. Of course everybody is welcome to post an answer!

Answer 1

MMA version: 12.0

I'll use the following images, found on the internet under open licence (if you know how to share the images themselves, please let me know). I originally had chosen two more but feature extraction did not work appropriately. I already rescaled them by chosing scaling factor such that the area of the triangle "left eye - right eye - mouth" was constant.

These are stored in the List imgs. Then, the features are extracted (in feats). Notice that I dropped "MouthInternalPoints" because there were of no real use but often contained duplicates which somehow complicates the process. The height and the width of the images are determined by the most restrictive image (here, the third one because the top of the head reaches the top of the frame).

feats = First@FacialFeatures[#, "Landmarks"] & /@ imgs;
feats = KeyDrop[#, "MouthInternalPoints"] & /@ feats;
xcenter = Mean[#["OutlinePoints"][[{1, -1}, 1]]] & /@ feats;
ycenter = #["NosePoints"][[1, 2]] & /@ feats;
{widthMax, heightMax} = Floor[Min /@ Transpose[
     2*(ImageDimensions /@ imgs - Transpose[{xcenter, ycenter}])]];
height = heightMax
width = Round[3/4*height]

Now that the image sizes are defined, the four vertices of the targeted frame are added and centered on the face outline (vertically) and the nose (horizontally):

addCorners[feat_, {width_, height_}] := 
 Block[{feattmp = feat, xcenter, ycenter, xmin, xmax, ymin, ymax},
  xcenter = feat["OutlinePoints"][[{1, -1}, 1]] // Mean;
  ycenter = feat["NosePoints"][[1, 2]];
  xmin = xcenter - width/2;
  xmax = xcenter + width/2;
  ymin = ycenter - height/2;
  ymax = ycenter + height/2;
  feattmp[
    "Corners"] = {{xmin, ymin}, {xmin, ymax}, {xmax, ymax}, {xmax, 
     ymin}};
  feattmp]

feats = addCorners[#, {width, height}] & /@ feats;

Table[HighlightImage[#, KeyValueMap[{Tooltip[#2, #1]} &, feats[[i]]]] &@imgs[[i]], 
   {i, Length@imgs}]

The images are then resized, the features points adapted accordingly:

 imgsCut = Table[ImageTrim[imgs[[i]], feats[[i]]["Corners"][[{1, 3}]]], {i, Length@imgs}]

points = Table[# - feat["Corners"][[1]] & /@ 
    Flatten[KeyValueMap[#2 &, feat], 1], {feat, feats}];
Table[HighlightImage[imgsCut[[i]], points[[i]]], {i, Length@imgsCut}]

Then, I reorder the image by "similarity" of the feature points (other criteria such as image distance might be more relevant):

g = WeightedAdjacencyGraph[DistanceMatrix[points]];
path = FindHamiltonianPath[g]
(* {5, 4, 2, 3, 1} *)

Then, the core:

(* i, starting image index; j, landing image index; step: smoothness of the animation *)
morph[i_, j_, step_] := 
 Block[{points1, points2, pointsMean, ldp, triangles1, triangles2, 
   mats, func, triangles, img1, img2, tab1, tab2},
  points1 = points[[i]];
  points2 = points[[j]];
  pointsMean = Mean /@ Transpose[{points[[i]], points[[j]]}];
  ldp = ListDensityPlot[ArrayPad[pointsMean, {0, {0, 1}}], 
    Mesh -> All, ColorFunction -> (White &)];
  trianglesOrder = First@Cases[ldp, Polygon[idx_] :> idx, Infinity];
  triangles1 = points1[[#]] & /@ trianglesOrder;
  triangles2 = points2[[#]] & /@ trianglesOrder;
  mats = Table[
    FindGeometricTransform[triangles2[[k]], triangles1[[k]], 
      TransformationClass -> "Affine"][[2, 1]], {k, 
     Length@triangles1}];
  func[{x_, y_}, coef_] := 
   Piecewise[
    Table[{TransformationFunction[
        coef*mats[[k]] + (1 - coef)*IdentityMatrix[3]][{x, y}], 
      inPolyQ[triangles1[[k]], {x, y}]}, {k, 1, 
      Length@triangles1}], {x, y}];
  triangles[coef_] := Map[func[#, 1 - coef] &, triangles1, {2}];
  img1[coef_] := 
   Graphics[{Texture[imgsCut[[i]]], 
     MapIndexed[
      Polygon[#, 
        VertexTextureCoordinates -> ({#[[1]]/width, #[[2]]/height} & /@
            triangles1[[First[#2]]])] &, triangles[coef]]}, 
    ImageSize -> {width, height}, PlotRangePadding -> None];
  img2[coef_] := 
   Graphics[{Texture[imgsCut[[j]]], 
     MapIndexed[
      Polygon[#, 
        VertexTextureCoordinates -> ({#[[1]]/width, #[[2]]/height} & /@
            triangles2[[First[#2]]])] &, triangles[1 - coef]]}, 
    ImageSize -> {width, height}, PlotRangePadding -> None];
  tab1 = Table[img1[coef], {coef, 0, 1, step}];
  tab2 = Table[img2[coef], {coef, 0, 1, step}];
  out = Table[
    Blend[{tab1[[-k]], tab2[[k]]}, (k - 1)/(Length[tab1] - 1)], {k, 1,
      Length[tab1]}];
  out]

inPolyQ[poly_, pt_] := Graphics`PolygonUtils`PointWindingNumber[poly, pt] =!= 0

Credits to:

• andre314 (https://mathematica.stackexchange.com/a/220452/18767) for an efficient piecewise affine transform
• rm -rf (https://mathematica.stackexchange.com/a/9417/18767) for inPolyQ

The principle is rather simple and common: perform a Delaunay mesh from the feature points, map each triangle from the origin image to the arrival image with an interpolation parameter (coef) that controls the smoothness of the transitions.

Then, the result:

pathCircular = path~Join~{path[[1]]};
out = ParallelTable[
   morph[pathCircular[[i]], pathCircular[[i + 1]], .1], {i, 
    Length[pathCircular] - 1}];

ListAnimate[Flatten@out]

(Poor quality due to 2Mb upload limitation)

Showing the mesh Simply add FaceForm[], EdgeForm[{Thick, Red}], Triangle /@ triangles[coef] in the Graphics of img1 and the same with 1-coef instead of coef in img2 to show the triangular mesh: