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:
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:
