*Valerie is a PhD student in computer science at Curtin University of Technology in Western Australia. Her main area of research is computer graphics, and her PhD topic is speech-driven facial animation. She can be reached at val@ marsh.cs.curtin.edu.au.*

To specify a morph, the animator defines a correspondence between the initial and final images. This is commonly done using points, triangles, or meshes. (Figure 1 shows the meshing process the Rmorf program uses.) Once the relationship has been defined, the textures in the images are blended from the initial image to the final image to produce the morph sequence.

The method of subdividing the images varies between different morphing programs. Triangulations can be specified by hand or generated by the morphing program. For example, in the morphing programs written by Mark Hall, Michael Kaufman, and Jay Windley, the user begins by specifying the triangles to be morphed in each image. Hall's program warps the triangular mesh of the input image to match the mesh of the output image. As the animation goes from the initial to the final image, the colors change from that of the first image to that of the second. Likewise, the user of Kaufman's program begins by defining common points in two images. For a human face, these points include the eyes, eyebrows, nose, and mouth. The program builds triangles from these points to create corresponding areas on the images. Windley's program does a linear blend of the contents of the triangles using barycentric geometry. One of the most automatic triangulation algorithms is the Delaunay triangulation, also known as the Voronoi tessellation. (See "Spatial Data and the Voronoi Tessellation," by Hrvoje Lukatela and John Russell, *DDJ*, December 1992.)

Morphing isn't always done with triangular patches. Douglas Smythe's program fits a Catmull-Rom spline through the x-coordinates of the control points (in the first pass) and the y-coordinates (second pass) to realize a piecewise continuous mapping function. Likewise, the Montage program for the Amiga (developed by Thomas Krehbiel and Kermit Woodall at Nova Design) uses cubic splines in a deformable mesh to define the input and output images. Cubic splines are better for defining curves in the image. Triangular patches have to be small and well chosen to give a similarly smooth result, but on the positive side, they are simpler to work with.

2-D morphing is a very manual process. It relies heavily on the experience and know-how of the animator. The choice of subject matter can make all the difference. Care should also be taken to use images of the subject captured from similar angles. Similarly, the position of parts of the subject should be closely matched.

- Is face connectivity maintained for all intermediate shapes?
- Are the intermediate shapes overly distorted?
- What restrictions exist on the original models?

Kent's system establishes a correspondence between the objects in three stages. First, each model is projected onto the surface of a unit sphere. Then the union of the two models' projections is taken using a modified version of the Weiler polygon-clipping algorithm. The merged topology is then projected back onto the original models. The results of Kent's approach are highly encouraging when working with convex and star-shaped objects. Future efforts will look at developing projections for general concave polyhedrons and polyhedrons with holes. They are also looking at giving the animator more control over the transformations. By using both topological and geometric information in their system, Kent maintains the integrity of the objects (in the Eulerian sense) while the output is intuitive in appearance. There's much less distortion of the shapes in this method.

Dave Bock, a visualization programmer at the North Carolina Supercomputing Center, takes a different approach to morphing by transforming objects using their mathematical descriptions. Bock generates volume data files that contain varying percentages of the two objects' data sets. The initial and final frames have 100 percent of one object and 0 percent of the other. He interpolates both volume data sets to define their respective objects at a common data threshold to use with an isosurface generator. Then a series of data files is created to hold the combined percentages of the initial data files at different stages of the morph. The number of data files depends on the number of in-between frames the user specifies. For each frame, volume data files are created containing increasing percentages of the final image and decreasing percentages of the initial object. Then the isosurface generator is used to create the geometric object for each frame, based on the volume data in the files. Bock is enhancing his system to accept geometrically defined objects as well as mathematically defined ones.

In all these 3-D morphing systems, the emphasis has been on getting the 3-D models defined--rendering and texture mapping is done later. This is one of the main differences between 2-D and 3-D morphing. Another difference is that 3-D morphing systems tend to be more automated than their 2-D counterparts.

Sorenson, P. "Morphing Magic." *Computer Graphics World* (January, 1992).

Wolberg, G. *Digital Image Warping*. Los Alamitos, CA: IEEE Computer Society Press, 1990.

A Shareware Morphing Program for MS-DOS

To this end, we're providing Rmorf, a shareware program for MS-DOS written by Richard Goedeken (see "Availability," page 5). Figure 1 shows the Rmorf user interface, which brings up two images, side by side, and allows you to specify corresponding mesh points in each image. The program then produces a series of in-between images, as many as the frame count you specify. You'll need a way to display the resulting image sequence, however. One possibility is to use Dave Mason's DTA or Autodesk's shareware TGAFLI. EXE program, both available on CompuServe, to chain together a sequence of Targa files into a Flic file that can then be displayed with a player such as the READFLIC program Jim Kent included with his article, "The Flic File Format" (*DDJ*, March 1993). In addition to Rmorf and the sample image files, we are providing TGAFLIC.EXEand READFLIC.EXE electronically (see "Availability," page 5).

Although industrial-strength morphing programs rely on interpolated splines and sophisticated antialiasing, Richard opts for a straightforward implementation and raw speed. Here's his description of the genesis of Rmorf:

Recently, I was browsing through *CompuServe Magazine* and noticed a small item about a shareware morphing program, so I logged on and downloaded it. The sparse documentation described the basics of image morphing, of which I was previously unaware. In running the program, it took me several attempts to produce a moderately decent morph.

After spending several hours with the software, I got to know its problems: It crashed frequently and was painfully slow.

I've done some work in computer graphics, so, over the following week, I considered the possibility of writing a better morphing program. I sketched out, in my mind and on paper, how two images could be morphed together. It wasn't as difficult as I had first imagined.

All I had to do was make a transfer mesh for each frame, which would be somewhere in between the first two meshes (depending upon the frame number), warp each image to the transfer mesh, and mix the colors in a given ratio (again depending upon the frame number). After approximately nine days of working in my spare time, I completed a rough copy of Rmorf. Two days later the finishing touches were done, and I released the program. The end result of these 11 days of intense work are about 3700 lines of source code--2400 in assembly language and the remaining 1300 in C.

Because the morphing portion is in assembler--and because it uses integer rather than floating-point math--Rmorf can calculate a frame in 7.2 seconds on my 33-MHz 386, instead of the 30--50 minutes required by other programs.

Rmorf currently supports only Targa files, but upcoming versions will include .FLI and .GIF support. The unregistered version of Rmorf supports 320x200 resolution; the registered version handles 1024x768.

Richard, a junior in high school, has written a number of other programs (a product-inventory system, a personal accounting package, and an SVGA game) which are available from his software company. He can be reached on CompuServe at 70304,1065.

--editors

Figure 1: Running the RMORF program. Note the mesh lines on the images which define a correspondence between one image and another.

The Canonical Implementation in C

The two-pass mesh warping algorithm implemented by Douglas Smythe at Industrial Light and Magic in 1989 has become the canonical implementation of morphing. In his classic book on digital image warping (*Digital Image Warping*, IEEE Computer Society Press, 1990), George Wolberg provides a lucid and thorough explanation of how the algorithm works, as well presenting his own version implemented in C. For this issue of *DDJ*, however, Wolberg has shared with us a new implementation that's more self-contained and optimized than the one in his book; see Listing One (page 92). For a complete description, refer to Wolberg's book and the documentation that accompanies the electronic version of this article; see "Availability," page 5. Here's a quick summary of Wolberg's description.

The algorithm relies on the fact that an arbitrary, one-pass, spatial transformation can be decomposed into a computationally cheaper two-pass operation. This approach stems from the seminal 1980 paper by Ed Catmull and Alvy Ray Smith, which is generally applicable to affine and perspective transformations on planar and nonplanar surfaces. Here, a 2-D resampling problem is replaced by two orthogonal 1-D resampling stages.

The user must specify two sets of control points. Both input and output images are thereby partitioned into a mesh of patches. Each patch delimits an image region over which a continuous mapping function applies. Mapping between both images now becomes a matter of transforming each patch onto its counterpart in the second image--known as *mesh warping*. Since the mapping function is defined only at these discrete points, it's necessary to determine the mapping function over all points to perform the warp. The patches can be fitted with a bivariate function to realize a piecewise continuous mapping function. Wolberg adds:

The benefit of using a mesh derives from the simplicity in interpolating the new positions of intermediate points (between the mesh points). A bilinear or bicubic function can be used. We use a Catmull-Rom cubic spline to implement bicubic interpolation here.

The algorithm requires that all four edges of each mesh be frozen. This means that the first and last rows and columns all remain intact throughout the warp. Wolberg continues:

The input includes a source image I1 and two meshes, M1 and M2. Mesh M1 is used to select landmark positions in I1, and M2 identifies their corresponding positions in the output image. In this manner, arbitrary points in I1 can be "pulled" to new positions. Although the use of a parametric mesh might seem to place unnecessary constraints on the positions of these points, a large class of useful transformations is possible. It is important, though, that the mesh not self-intersect in order to avoid the image from folding upon itself.

The algorithm's first pass is responsible for resampling each row independently. It maps all (u,v) points in the source image I1 to their (x,v) coordinates in the intermediate image, thereby positioning each input point into its proper output column. The intermediate image is that whose x-coordinates are the same as those in I2 and whose y-coordinates are taken from I1. The second pass then resamples each column in the intermediate image, mapping every (x,v) to its final (x,y) position. Each point now lies in its proper row as well as its column.

The implementation presented here is for gray-scale images only. It is straightforward to extend the program to handle three-channel images (like RGB color images), by handling each channel separately.

Wolberg's implementation is in the classic UNIX style of a rudimentary command-line user interface. Wolberg explains:

The code is missing a program to help the user create and edit meshes interactively. A good mesh editor is a critical component to any mesh warping program. Such code falls outside of the scope of this presentation. A more full-featured implementation would allow the user to control the cross-dissolve schedule at each mesh point, as well as its position. This permits the intensities in different regions of the image to interpolate at different rates.

--editors

_MORPHING IN 2D AND 3D_ by Valerie Hall Source code written by George Wolberg and accompanies the sidebar entitled "The Canonical Implementation in C" /* ======================================================================== * meshwarp.h -- Header file for meshwarp.c. (C) 1993 by George Wolberg. * ======================================================================*/ #include <stdio.h> #define BW 0 #define MESH 1 #define MAX(A,B) ((A) > (B) ? (A) : (B)) #define MIN(A,B) ((A) < (B) ? (A) : (B)) typedef unsigned char uchar; typedef struct { /* image data structure */ int width; /* image width (# cols) */ int height; /* image height (# rows) */ void *ch[2]; /* pointers to channels */ } imageS, *imageP; extern void meshWarp(); /* extern decls for funcs in meshwarp.c */ extern void resample(); extern imageP readImage(); /* extern decls for funcs in util.c */ extern int saveImage(); extern imageP allocImage(); extern void freeImage(); /* ======================================================================= * morph.c - Generate a metamorphosis sequence. (C) 1993 by George Wolberg. * =======================================================================*/ #include "meshwarp.h" /*------- main: Collect user parameters and pass them to morph() ------- */ main(argc, argv) int argc; char **argv; { int nframes; char name[10]; imageP I1, I2; imageP M1, M2; if(argc != 7) /* make sure user invokes this program properly */ { fprintf(stderr, "Usage: morph src.bw dst.bw src.XY dst.XY frames name\n"); exit(1); } /*----------- read input image and meshes --------------*/ I1 = readImage(argv[1], BW); /* source image */ I2 = readImage(argv[2], BW); /* target image */ M1 = readImage(argv[3], MESH); /* source mesh */ M2 = readImage(argv[4], MESH); /* target mesh */ nframes = atoi(argv[5]); /* # frames */ strcpy(name, argv[6]); /* out basename */ /*----------- call morph -------------------------------*/ morph(I1, I2, M1, M2, nframes, name); } /* ------------------------------------------------------------------------ * morph: Generate a morph sequence of frames between images I1 and I2. * Correspondence points among I1 and I2 are given in meshes M1 and M2. * nframes frames are generated (including I1 and I2). The output is stored * in files "basename_xxx.bw" where xxx are sequential 3-digit frame numbers. *------------------------------------------------------------------------*/ void morph(I1, I2, M1, M2, nframes, basename) imageP I1, I2, M1, M2; int nframes; char * basename; { int i, j, totalI, totalM; double w1, w2; char name[20]; uchar *p1, *p2, *p3; float *x1, *y1, *x2, *y2, *x3, *y3; imageP I3, Iw1, Iw2, M3; /* allocate space for tmp images and mesh */ M3 = allocImage(M1->width, M1->height, MESH); I3 = allocImage(I1->width, I1->height, BW); Iw1 = allocImage(I1->width, I1->height, BW); Iw2 = allocImage(I1->width, I1->height, BW); /* eval total number of points in mesh (totalM) and image (totalI) */ totalM = M1->width * M1->height; totalI = I1->width * I1->height; /* copy 1st frame to basename_000.bw*/ sprintf(name, "%s_000.bw", basename); saveImage(I1, name, BW); printf("Finished Frame 0\n"); for(i=1; i<nframes-1; i++) { /* M3 <- linearly interpolate between M1 and M2 */ w2 = (double) i / (nframes-1); w1 = 1. - w2; /* linearly interpolate M3 grid */ x1 = (float *) M1->ch[0]; y1 = (float *) M1->ch[1]; x2 = (float *) M2->ch[0]; y2 = (float *) M2->ch[1]; x3 = (float *) M3->ch[0]; y3 = (float *) M3->ch[1]; for(j=0; j<totalM; j++) { x3[j] = x1[j]*w1 + x2[j]*w2; y3[j] = y1[j]*w1 + y2[j]*w2; } /* warp I1 and I2 according to grid M3 */ meshWarp(I1, M1, M3, Iw1); meshWarp(I2, M2, M3, Iw2); /* cross-dissolve warped images Iw1 and Iw2 */ p1 = (uchar *) Iw1->ch[0]; p2 = (uchar *) Iw2->ch[0]; p3 = (uchar *) I3->ch[0]; for(j=0; j<totalI; j++) p3[j] = p1[j]*w1 + p2[j]*w2; /* save frame into file */ sprintf(name, "%s_%03d.bw", basename, i); saveImage(I3, name, BW); printf("Finished Frame %d\n", i); } /* copy last frame to basename_xxx.bw */ sprintf(name, "%s_%03d.bw", basename, i); saveImage(I2, name, BW); printf("Finished Frame %d\n", i); } /* =========================================================================== * catmullRom.c - Catmull-Rom interpolating spline. (C) 1993 by George Wolberg * =======================================================================*/ #include "meshwarp.h" /* ------------------------------------------------------------------------ * catmullRom: Compute a Catmull-Rom spline passing thru the len1 points in * arrays x1, y1, where y1 = f(x1). len2 positions on the spline are to be * Their positions are given in x2. The spline values are stored in y2. *------------------------------------------------------------------------*/ void catmullRom(x1, y1, len1, x2, y2, len2) float *x1, *y1, *x2, *y2; int len1, len2; { int i, j, dir, j1, j2; double x, dx1, dx2; double dx, dy, yd1, yd2, p1, p2, p3; double a0y, a1y, a2y, a3y; /* find direction of monotonic x1; skip ends */ if(x1[0] < x1[1]) /* increasing */ { if(x2[0]<x1[0] || x2[len2-1]>x1[len1-1]) dir=0; else dir = 1; } else /* decreasing */ { if(x2[0]>x1[0] || x2[len2-1]<x1[len1-1]) dir=0; else dir = -1; } if(dir == 0) /* error */ { printf("catmullRom: Output x-coord out of range of input\n"); return; } /* p1 is first endpoint of interval * p2 is resampling position * p3 is second endpoint of interval * j is input index for current interval */ if(dir==1) p3 = x2[0] - 1; /* force coefficient initialization */ else p3 = x2[0] + 1; for(i=0; i<len2; i++) { p2 = x2[i]; /* check if in new interval */ if( (dir== 1 && p2>p3 ) || (dir== -1 && p2<p3 )) { if(dir) /* find the interval which contains p2 */ { for(j=0; j<len1 && p2>x1[j]; j++) ; if(p2 < x1[j]) j--; } else { for(j=0; j<len1 && p2<x1[j]; j++) ; if(p2 > x1[j]) j--; } p1 = x1[j]; /* update 1st endpt */ p3 = x1[j+1]; /* update 2nd endpt */ /* clamp indices for endpoint interpolation */ j1 = MAX(j-1, 0); j2 = MIN(j+2, len1-1); /* compute spline coefficients */ dx = 1.0 / (p3 - p1); dx1 = 1.0 / (p3 - x1[j1]); dx2 = 1.0 / (x1[j2] - p1); dy = (y1[j+1] - y1[ j ]) * dx; yd1 = (y1[j+1] - y1[ j1]) * dx1; yd2 = (y1[j2 ] - y1[ j ]) * dx2; a0y = y1[j]; a1y = yd1; a2y = dx * ( 3*dy - 2*yd1 - yd2); a3y = dx*dx*(-2*dy + yd1 + yd2); } /* use Horner's rule to calculate cubic polynomial */ x = p2 - p1; y2[i] = ((a3y*x + a2y)*x + a1y)*x + a0y; } } /* ======================================================================== * meshwarp.c -- Mesh warping program. Copyright (C) 1993 by George Wolberg. * =======================================================================*/ #include "meshwarp.h" /* ------------------------------------------------------------------------- * meshWarp: Warp I1 with correspondence points given in meshes M1 and M2. * Result goes in I2. Based on Douglas Smythe's algorithm at ILM. *------------------------------------------------------------------------*/ void meshWarp(I1, M1, M2, I2) imageP I1, I2, M1, M2; { int I_w, I_h, M_w, M_h; int x, y, u, v, n; float *x1, *y1, *x2, *y2; float *xrow, *yrow, *xcol, *ycol, *coll, *indx, *map; uchar *src, *dst; imageP Mx, My, I3; I_w = I1->width; I_h = I1->height; M_w = M1->width; M_h = M1->height; /* alloc enough memory for a scanline along the longest dimension */ n = MAX(I_w, I_h); indx = (float *) malloc(n * sizeof(float)); /* should check if err */ xrow = (float *) malloc(n * sizeof(float)); yrow = (float *) malloc(n * sizeof(float)); map = (float *) malloc(n * sizeof(float)); /* create table of x-intercepts for source mesh's vert splines */ Mx = allocImage(M_w, I_h, MESH); for(y=0; y < I_h; y++) indx[y] = y; for(u=0; u < M_w; u++) /* visit each vert spline */ { /* store col as row for spline function */ xcol = (float *) M1->ch[0] + u; ycol = (float *) M1->ch[1] + u; coll = (float *) Mx->ch[0] + u; /* scan-convert vert splines */ for(v=0; v < M_h; v++, xcol+=M_w) xrow[v] = *xcol; for(v=0; v < M_h; v++, ycol+=M_w) yrow[v] = *ycol; catmullRom(yrow, xrow, M_h, indx, map, I_h); /* store resampled row back into column */ for(y=0; y < I_h; y++, coll+=M_w) *coll = map[y]; } /* create table of x-intercepts for dst mesh's vert splines */ for(u=0; u < M_w; u++) /* visit each vert spline */ { /* store column as row for spline fct */ xcol = (float *) M2->ch[0] + u; ycol = (float *) M2->ch[1] + u; coll = (float *) Mx->ch[1] + u; /* scan-convert vert splines */ for(v=0; v < M_h; v++, xcol+=M_w) xrow[v] = *xcol; for(v=0; v < M_h; v++, ycol+=M_w) yrow[v] = *ycol; catmullRom(yrow, xrow, M_h, indx, map, I_h); /* store resampled row back into column */ for(y=0; y < I_h; y++, coll+=M_w) *coll = map[y]; } /*------------ first pass: warp x using tables in Mx --------*/ I3 = allocImage(I_w, I_h, BW); x1 = (float *) Mx->ch[0]; x2 = (float *) Mx->ch[1]; src = (uchar *) I1->ch[0]; dst = (uchar *) I3->ch[0]; for(x=0; x < I_w; x++) indx[x] = x; for(y=0; y < I_h; y++) { /* fit spline to x-intercepts; resample over all cols */ catmullRom(x1, x2, M_w, indx, map, I_w); /* resample source row based on map */ resample(src, I_w, 1, map, dst); /* advance pointers to next row */ src += I_w; dst += I_w; x1 += M_w; x2 += M_w; } freeImage(Mx); /* create table of y-intercepts for intermediate mesh's hor splines */ My = allocImage(I_w, M_h, MESH); x1 = (float *) M2->ch[0]; y1 = (float *) M1->ch[1]; y2 = (float *) My->ch[0]; for(x=0; x < I_w; x++) indx[x] = x; for(v=0; v < M_h; v++) /* visit each horz spline */ { /* scan-convert horz splines */ catmullRom(x1, y1, M_w, indx, y2, I_w); x1 += M_w; /* advance pointers to next row */ y1 += M_w; y2 += I_w; } /* create table of y-intercepts for dst mesh's horz splines */ x1 = (float *) M2->ch[0]; y1 = (float *) M2->ch[1]; y2 = (float *) My->ch[1]; for(v=0; v < M_h; v++) /* visit each horz spline */ { /* scan-convert horz splines */ catmullRom(x1, y1, M_w, indx, y2, I_w); x1 += M_w; /* advance pointers to next row */ y1 += M_w; y2 += I_w; } /*----------------------- second pass: warp y ------------------*/ src = (uchar *) I3->ch[0]; dst = (uchar *) I2->ch[0]; for(y=0; y < I_h; y++) indx[y] = y; for(x=0; x < I_w; x++) { /* store column as row for spline fct */ xcol = (float *) My->ch[0] + x; ycol = (float *) My->ch[1] + x; for(v=0; v < M_h; v++, xcol+=I_w) xrow[v] = *xcol; for(v=0; v < M_h; v++, ycol+=I_w) yrow[v] = *ycol; /* fit spline to y-intercepts; resample over all rows */ catmullRom(xrow, yrow, M_h, indx, map, I_h); /* resample source column based on map */ resample(src, I_h, I_w, map, dst); /* advance pointers to next column */ src++; dst++; } freeImage(My); freeImage(I3); free((char *) indx); free((char *) xrow); free((char *) yrow); free((char *) map); } /* ------------------------------------------------------------------------ * resample: Resample the len elements of src (with stride offst) into dst * according to the spatial mapping given in xmap. Perform linear interpola- * tion for magnification and box filtering (unweighted averaging) for * minification. Based on Fant's algorithm (IEEE Comp. Graphics & Appl. 1/86) *------------------------------------------------------------------------*/ void resample(src, len, offst, xmap, dst) uchar *src, *dst; float *xmap; int len, offst; { int u, x, v0, v1; double val, sizfac, inseg, outseg, acc, inpos[1024]; /* precompute input index for each output pixel */ for(u=x=0; x<len; x++) { while(xmap[u+1]<x) u++; inpos[x] = u + (double) (x-xmap[u]) / (xmap[u+1]-xmap[u]); } inseg = 1.0; outseg = inpos[1]; sizfac = outseg; acc = 0.; v0 = *src; src += offst; v1 = *src; src += offst; for(u=1; u<len; ) { val = inseg*v0 + (1-inseg)*v1; if(inseg < outseg) { acc += (val * inseg); outseg -= inseg; inseg = 1.0; v0 = v1; v1 = *src; src += offst; } else { acc += (val * outseg); acc /= sizfac; *dst = (int) MIN(acc, 0xff); dst += offst; acc = 0.; inseg -= outseg; outseg = inpos[u+1] - inpos[u]; sizfac = outseg; u++; } } }