_MORPHING IN 2D AND 3D_
by Valerie Hall
Source code written by George Wolberg and accompanies the 
sidebar entitled "The Canonical Implementation in C"

[LISTING ONE]

/* ==========================================================================
 * 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++;
                }
        }
}