/*
 * Copyright (c) 2002 Mark Salyzyn
 * All rights reserved.
 *
 * TERMS AND CONDITIONS OF USE
 *
 * Redistribution and use in source form, with or without modification, are
 * permitted provided that redistributions of source code must retain the
 * above copyright notice, this list of conditions and the following
 * disclaimer.
 *
 * This software is provided `as is' by Mark Salyzyn and any express or implied
 * warranties, including, but not limited to, the implied warranties of
 * merchantability and fitness for a particular purpose, are disclaimed. In no
 * event shall Mark Salyzyn be liable for any direct, indirect, incidental,
 * special, exemplary or consequential damages (including, but not limited to,
 * procurement of substitute goods or services; loss of use, data, or profits;
 * or business interruptions) however caused and on any theory of liability,
 * whether in contract, strict liability, or tort (including negligence or
 * otherwise) arising in any way out of the use of this driver software, even
 * if advised of the possibility of such damage.
 *
 *      Name: pnmstitch.c
 *  Description: Automated panoramic stitcher, with the further
 *  constraint that the right hand side of the right hand image is `fixed'
 *  per-se, after offset evaluation and only the left hand side of the
 *  right hand image is mangled.  Thus, the algorithm is optimized for
 *  stitching a right hand image to the right hand half (half being a lose
 *  term) of the left hand image.
 *  Author: Mark Salyzyn June 2002
 *
 *  ToDo:
 *      - Add and refine the videorbits algorithm
 *      - speedups in matching algorithm
 *      - refinement in matching algorithm
 *      - Add RotateSliver transform algorithm (reuse SliverMatch,
 *        and rotate the right image).
 */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
/*
#include <stdint.h>
*/
#define INTPTR_MIN        (-0x7fffffff - 1)
#define INTPTR_MAX        0x7fffffff
#define UINTPTR_MAX       0xffffffff

#include "pam.h"
#include "shhopt.h"

/*
 *  Defines
 */
#if (!defined(UNREFERENCED_PARAMETER))
# if (defined(__GNUC__))
#  define UNREFERENCED_PARAMETER(x)
# elif (defined(__USLC__) || defined(_M_XENIX))
#  define UNREFERENCED_PARAMETER(x) ((x)=(x))
# else
#  define UNREFERENCED_PARAMETER(x) (x)
# endif
#endif

#if (!defined(__cplusplus))
# if (!defined(false))
#  define false ((bool)0)
# endif
# if (!defined(true))
#  define true ((bool)1)
# endif
#endif

/*
 *  Structures
 */

/*
 *  Image structure
 */
typedef struct {
    char     * name;
    struct pam pam;
    tuple   ** tuple;
} Image;

/*
 *  Output class
 */
typedef struct output {
    /* name */
    char    *Name;
    /* methods */
    bool   (*Alloc)(struct output * me, char * file, unsigned width, unsigned height, struct pam * prototype);
    void   (*DeAlloc)(struct output * me);
    tuple *(*Row)(struct output * me, unsigned row);
    void   (*FlushRow)(struct output * me, unsigned row);
    void   (*FlushImage)(struct output * me);
    /* data */
    Image   *image;
    void    *extra;
} Output;

extern Output OutputMethods[];

/*
 *  Stitching class
 */
typedef struct stitcher {
    /* name */
    char  *Name;
    /* methods */
    bool (*Alloc)(struct stitcher *me);
    void (*DeAlloc)(struct stitcher *me);
    /* Inform transformation parameter routines about constraints */
    void (*Constrain)(struct stitcher *me, int x, int y, int width, int height);
    /* transformation parameter determination */
    bool (*Match)(struct stitcher *me, Image * Left, Image * Right);
    /* transformation methods */
    int  (*XLeft)(struct stitcher *me, int x, int y);
    int  (*YLeft)(struct stitcher *me, int x, int y);
    int  (*XRight)(struct stitcher *me, int x, int y);
    int  (*YRight)(struct stitcher *me, int x, int y);
    /* Output methods */
    void (*Output)(struct stitcher *me, FILE * fp);
    /* data */
    int    x, y, width, height;
    float *parms;
} Stitcher;

extern Stitcher StitcherMethods[];

/*
 *  Prototypes
 */
static int pnmstitch(char * left, char * right, char * out,
                     int x, int y, int width, int height,
                     char * stitcher, char * filter);

struct cmdlineInfo {
    /* All the information the user supplied in the command line,
       in a form easy for the program to use.
    */
    char *leftFilespec;  /* '-' if stdin */
    char *rightFilespec;
    char *outputFilespec;
    char *stitcher;
    char *filter;
    int width;
    int height;
    int xrightpos;
    int yrightpos;
    unsigned int verbose;
};




static void
parseCommandLine ( int argc, char ** argv,
                   struct cmdlineInfo *cmdlineP ) {
/*----------------------------------------------------------------------------
   parse program command line described in Unix standard form by argc
   and argv.  Return the information in the options as *cmdlineP.  

   If command line is internally inconsistent (invalid options, etc.),
   issue error message to stderr and abort program.

   Note that the strings we return are stored in the storage that
   was passed to us as the argv array.  We also trash *argv.
-----------------------------------------------------------------------------*/
    optEntry *option_def = malloc( 100*sizeof( optEntry ) );
        /* Instructions to optParseOptions3 on how to parse our options.
         */
    optStruct3 opt;

    unsigned int option_def_index;

    char *outputOpt;
    unsigned int widthSpec, heightSpec, outputSpec, 
        xrightposSpec, yrightposSpec, stitcherSpec, filterSpec;

    option_def_index = 0;   /* incremented by OPTENT3 */
    OPTENT3(0, "width",       OPT_UINT,   &cmdlineP->width, 
            &widthSpec,         0);
    OPTENT3(0, "height",      OPT_UINT,   &cmdlineP->height, 
            &heightSpec,        0);
    OPTENT3(0, "verbose",     OPT_FLAG,   NULL,                  
            &cmdlineP->verbose, 0 );
    OPTENT3(0, "output",      OPT_STRING, &outputOpt, 
            &outputSpec,        0);
    OPTENT3(0, "xrightpos",   OPT_UINT,   &cmdlineP->xrightpos, 
            &xrightposSpec,     0);
    OPTENT3(0, "yrightpos",   OPT_UINT,   &cmdlineP->yrightpos, 
            &yrightposSpec,     0);
    OPTENT3(0, "stitcher",    OPT_STRING, &cmdlineP->stitcher, 
            &stitcherSpec,      0);
    OPTENT3(0, "filter",      OPT_STRING, &cmdlineP->filter, 
            &filterSpec,        0);

    opt.opt_table = option_def;
    opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
    opt.allowNegNum = FALSE;  /* We have no parms that are negative numbers */

    optParseOptions3( &argc, argv, opt, sizeof(opt), 0);
        /* Uses and sets argc, argv, and some of *cmdlineP and others. */

    if (!widthSpec)
        cmdlineP->width = INTPTR_MAX;
    if (!heightSpec)
        cmdlineP->height = INTPTR_MAX;
    if (!xrightposSpec)
        cmdlineP->xrightpos = INTPTR_MAX;
    if (!yrightposSpec)
        cmdlineP->yrightpos = INTPTR_MAX;
    if (!stitcherSpec)
        cmdlineP->stitcher = "BiLinearSliver";
    if (!filterSpec)
        cmdlineP->filter = "LineAtATime";


    if (argc-1 > 3)
        pm_error("Program takes at most three arguments: left, right, and "
                 "output file specifications.  You specified %d", argc-1);
    else {
        if (argc-1 == 0) {
            cmdlineP->leftFilespec = "-";
            cmdlineP->rightFilespec = "-";
        } else if (argc-1 == 1) {
            cmdlineP->leftFilespec = "-";
            cmdlineP->rightFilespec = argv[1];
        } else {
            cmdlineP->leftFilespec = argv[1];
            cmdlineP->rightFilespec = argv[2];
        }
        if (argc-1 == 3 && outputSpec)
            pm_error("You cannot specify --output and also name the "
                     "output file with the 3rd argument.");
        else if (argc-1 == 3)
            cmdlineP->outputFilespec = argv[3];
        else if (outputSpec) 
            cmdlineP->outputFilespec = outputOpt;
        else
            cmdlineP->outputFilespec = "-";
    }
}




static char minus[] = "-";
static int  verbose;

/*
 *  Parse the command line, call pnmstitch to perform work.
 */
int
main (int argc, char **argv)
{
    struct cmdlineInfo cmdline;

    pnm_init(&argc, argv);

    parseCommandLine(argc, argv, &cmdline);

    verbose = cmdline.verbose;

    return pnmstitch (cmdline.leftFilespec, 
                      cmdline.rightFilespec, 
                      cmdline.outputFilespec, 
                      cmdline.xrightpos, 
                      cmdline.yrightpos, 
                      cmdline.width, 
                      cmdline.height, 
                      cmdline.stitcher, 
                      cmdline.filter);
} /* main() - end */

/*
 *  allocate a clear image structure.
 */
static Image *
allocate_image(void)
{
    Image * retVal = malloc((unsigned)sizeof(Image));

    if (retVal == NULL) {
        return NULL;
    }
    memset (retVal, 0, (unsigned)sizeof(Image));
    return retVal;
} /* allocate_image() - end */

/*
 *  free a clear image structure.
 */
static void
free_image(Image * image)
{
    if (image->name) {
        free (image->name);
    }
    if (image->tuple) {
        pnm_freepamarray(image->tuple, &image->pam);
    }
    free (image);
} /* free_image() - end */

/*
 *  Initialize a tuple array with the contents of the image file.
 */
static Image *
readinit(char * name)
{
    FILE * ifP;
    Image * image = allocate_image();

    if (image == NULL) {
        return NULL;
    }
    if (strcmp(name,minus) == 0) {
        ifP = stdin;
        image->name = strdup("<stdin>");
    } else {
        /* list of possible extensions to read from */
        char  * list[] = { ".pnm", ".pam", ".pgm", ".pbm", ".ppm" };
        char ** entry = list - 1;

        image->name = strdup(name);
        for (;;) {
            ifP = fopen(image->name,"rb");
            if (ifP != NULL) {
                break;
            }
            free (image->name);
            image->name = NULL;
            if ((++entry >= &list[sizeof(list)/sizeof(list[0])])
                || ((image->name = malloc(strlen(name)+strlen(*entry)+1)) == 
                    NULL)) {
                break;
            }
            strcat(strcpy(image->name,name),*entry);
        }
        if (ifP == NULL) {
            pm_error ("Failed to open \"%s\"\n", name);
        }
    }
    image->tuple = 
        pnm_readpam(ifP, &(image->pam), (unsigned)sizeof(image->pam));
    fclose (ifP);

    return image;
} /* readinit() - end */

/*
 *  Prepare an image to be output.
 */
static bool
writeinit(Image * image)
{
    image->pam.file = pm_openw(image->name);

    pnm_writepaminit(&image->pam);

    if (strcmp(image->name,minus) == 0) {
        free (image->name);
        image->name = strdup("<stdout>");
    } 

    return true;
} /* writeinit() - end */




/*
   Return the magnitude of the difference between the specified region
   in the left image and the same region in the y image.  The magnitude
   is defined as the sum of the squares of the differences between 
   intensities of each of the 3 colors over all the pixels.
*/
static unsigned long
sub_tuple(Image * left, int lx, int ly, Image * right, int rx, int ry, 
          int width, int height)
{
    unsigned long total = 0;
    tuple      ** lROW = &(left->tuple)[ly];
    tuple      ** rROW = &(right->tuple)[ry];
    unsigned      row;

    for (row = 0; row < height; ++row) {
        tuple * lCOLUMN = &(*(lROW++))[lx];
        tuple * rCOLUMN = &(*(rROW++))[rx];
        unsigned column;
        for (column = 0; column < width; ++column) {
            tuple lPLANE = &(*(lCOLUMN++))[0];
            tuple rPLANE = &(*(rCOLUMN++))[0];
            unsigned plane;
            for (plane = 0; plane < left->pam.depth; ++plane) {
                long difference = *(lPLANE++) - *(rPLANE++);
                total += (unsigned long)(difference * difference);
            }
        }
    }
    return total;
} /* sub_tuple() - end */

/*
 *  Generate a recent sorted histogram of best matches.
 */
typedef struct {
    unsigned long total;
    int           x;
    int           y;
} Best;
/* Arbitrary, except 9 points surrounding a hot spot plus one layer more */
#define NUM_BEST 9

static Best *
allocate_best(void)
{
    Best * retVal = malloc((unsigned)sizeof(Best)*NUM_BEST);
    int    i;

    if (retVal == NULL) {
        return NULL;
    }
    for (i = 0; i < NUM_BEST; ++i) {
        retVal[i].total = UINTPTR_MAX;
        retVal[i].x = INTPTR_MAX;
        retVal[i].y = INTPTR_MAX;
    }
    return retVal;
} /* allocate_best() - end */

static void
update_best(Best * best, unsigned long total, int x, int y)
{
    int i;
    if (best[NUM_BEST-1].total <= total) {
        return;
    }
    for (i = NUM_BEST - 1; i > 0; --i) {
        if (best[i-1].total < total) {
            break;
        }
        best[i] = best[i-1];
    }
    best[i].total = total;
    best[i].x = x;
    best[i].y = y;
}

static void
pr_best(Best * best)
{
    int i;
    fprintf (stderr,"%p", best);
    for (i = 0; i < NUM_BEST; ++i) {
        fprintf (stderr, " (%d,%d)%lu", best[i].x, best[i].y, best[i].total);
    }
} /* pr_best() - end */



/*
 *  The wrapper for both the Output and the Stitcher algorithms.
 */

static void
determineMaskCorners(Stitcher * Stitch, 
                     Image * Left, Image * Right,
                     int xp[], int yp[]) {

    /* Determine the mask corners for both Left and Right images */

    int i;

    xp[0] = xp[1] = xp[4] = xp[5] = 0;
    yp[0] = yp[2] = yp[4] = yp[6] = 0;
    xp[2] = xp[3] = Left->pam.width;
    yp[1] = yp[3] = Left->pam.height;
    xp[6] = xp[7] = Right->pam.width;
    yp[5] = yp[7] = Right->pam.height;
    for (i = 0; i < 8; ++i) {
        int   xx, yy, count = 65536;
        int x,y;
        int (*X)(Stitcher *me, int x, int y);
        int (*Y)(Stitcher *me, int x, int y);

        if (i < 4) {
            X = Stitch->XLeft;
            Y = Stitch->YLeft;
        } else {
            X = Stitch->XRight;
            Y = Stitch->YRight;
        }
        x = xp[i];
        y = yp[i];
        /* will not work if rotated 90o or if gain > 10 */
        do {
            xx = (*X)(Stitch, xp[i], yp[i]) - x;
            if (xx < 0) {
                if (xx > -100) {
                    ++xp[i];
                } else {
                    xp[i] -= xx / 10;
                }
            } else if (xx > 0) {
                if (xx < 100) {
                    --xp[i];
                } else {
                    xp[i] -= xx / 10;
                }
            }
            yy = (*Y)(Stitch, xp[i], yp[i]) - y;
            if (yy < 0) {
                if (yy > -100) {
                    ++yp[i];
                } else {
                    yp[i] -= yy / 10;
                }
            } else if (yy > 0) {
                if (yy < 100) {
                    --yp[i];
                } else {
                    yp[i] -= yy / 10;
                }
            }
        } while (((xx != 0) || (yy != 0)) && (--count != 0));
    }
    if (verbose) {
        (*(Stitch->Output))(Stitch, stderr);
        if (verbose > 2) {
            static char quotes[] = "'\0\0\"\0\0'\"\0\"\"";
            fprintf (stderr, " Left:");
            for (i = 0; i < 8; ++i) {
                if (i == 4) {
                    fprintf (stderr, "\n Right:");
                }
                fprintf (stderr, " x%s,y%s=%d,%d",
                  &quotes[(i%4)*3], &quotes[(i%4)*3],
                  xp[i], yp[i]);
            }
        }
    }
}



static void
calculateXyWidthHeight(int xp[], int yp[], 
                       int * const xP, int * const yP,
                       int * const widthP, int * const heightP) {
    
    /* Calculate generic x,y left top corner, and the width and height */

    int x, y;
    int width, height;
    int i;

    x = xp[0];
    y = yp[0];
    width = height = 0;
    for (i = 1; i < 8; ++i) {
        if (xp[i] < x) {
            width += x - xp[i];
            x = xp[i];
        } else if ((x + width) < xp[i]) {
            width = xp[i] - x;
        }
        if (yp[i] < y) {
            height += y - yp[i];
            y = yp[i];
        } else if ((y + height) < yp[i]) {
            height = yp[i] - y;
        }
    }
    *xP = x; *yP = y;
    *widthP = width; *heightP = height;
}



static void
printPlan(int xp[], int yp[], Image * Left, Image * Right) {
    /* Calculate Left image transformed bounds */
    int X, Y, W, H;
    int i;

    X = xp[0];
    Y = yp[0];
    W = H = 0;
    for (i = 1; i < 4; ++i) {
        if (xp[i] < X) {
            W += X - xp[i];
            X = xp[i];
        } else if ((X + W) < xp[i]) {
            W = xp[i] - X;
        }
        if (yp[i] < Y) {
            H += Y - yp[i];
            Y = yp[i];
        } else if ((Y + H) < yp[i]) {
            H = yp[i] - Y;
        }
    }
    fprintf (stderr,
             "%s[%u,%u=>%d,%d](%d,%d)",
             Left->name, Left->pam.width, Left->pam.height,
             W, H, X, Y);
    X = xp[i];
    Y = yp[i];
    W = H = 0;
    for (++i; i < 8; ++i) {
        if (xp[i] < X) {
            W += X - xp[i];
            X = xp[i];
        } else if ((X + W) < xp[i]) {
            W = xp[i] - X;
        }
        if (yp[i] < Y) {
            H += Y - yp[i];
            Y = yp[i];
        } else if ((Y + H) < yp[i]) {
            H = yp[i] - Y;
        }
    }
    fprintf (stderr,
             "+%s[%u,%u=>%d,%d](%d,%d)",
             Right->name, Right->pam.width, Right->pam.height,
             W, H, X, Y);
}


static int
pnmstitch(char * const left, char * const right, char * const out, 
          int const reqx, int const reqy, 
          int const reqWidth, int const reqHeight, char * const stitcher,
          char * const filter)
{
    Image    * Left;
    Image    * Right;
    Output   * Out = &OutputMethods[0];
    Stitcher * Stitch = &StitcherMethods[0];
    unsigned   row;
    int        xp[8], yp[8];
    int        width, height;
    int        x, y;

    /* Select output filter */
    for (; Out->Name != (char *)NULL; ++Out) {
        if (strcasecmp (Out->Name, filter) == 0) {
            break;
        }
    }
    if (Out->Name == (char *)NULL) {
        fprintf (stderr,
          "Unknown output filter %s. available filters:\n", filter);
        for (Out = &OutputMethods[0]; Out->Name != (char *)NULL; ++Out) {
            fprintf (stderr, "\t%s\n", Out->Name);
        }
        return -1;
    }
    if (verbose) {
        fprintf (stderr, "Selected %s output filter\n", Out->Name);
    }
    /* Select stitcher type */
    for (; Stitch->Name != (char *)NULL; ++Stitch) {
        if (strcasecmp (Stitch->Name, stitcher) == 0) {
            break;
        }
    }
    if (Stitch->Name == (char *)NULL) {
        fprintf (stderr,
          "Unknown stitcher algorithm %s. Available stitchers:\n",
          stitcher);
        for (Stitch = &StitcherMethods[0]; 
             Stitch->Name != (char *)NULL; 
             ++Stitch) {
            fprintf (stderr, "\t%s\n", Stitch->Name);
        }
        return -2;
    }
    if (verbose) {
        fprintf (stderr, "Selected %s stitcher algorithm\n",
          Stitch->Name);
    }
    /* Left hand image read into memory */
    if ((Left = readinit(left)) == NULL) {
        return -3;
    }
    /* Right hand image read into memory */
    if ((Right = readinit(right)) == NULL) {
        return -4;
    }
    if ((*(Stitch->Alloc))(Stitch) == false) {
        return -5;
    }
    (*(Stitch->Constrain))(Stitch, reqx, reqy, reqWidth, reqHeight);
    if ((*(Stitch->Match))(Stitch, Left, Right) == false) {
        return -6;
    }

    determineMaskCorners(Stitch, Left, Right, xp, yp);

    /* Output the combined images */

    calculateXyWidthHeight(xp, yp, &x, &y, &width, &height);

    if (verbose) 
        printPlan(xp, yp, Left, Right);

    if ((*(Out->Alloc))(Out, out, width, height, &Left->pam) == false) {
        return -7;
    }
    if (verbose) {
        fprintf (stderr,
          "=%s[%u,%u=>%d,%d](%d,%d)\n",
          Out->image->name, Out->image->pam.width,
          Out->image->pam.height, width, height, x, y);
    }
    /*
     *  We scale the overlap of the left and right images, we need to
     * discover and hold on to the left edge of the right image to
     * determine the rate at which we blend. Most (7/8) of the blending
     * occurs in the first half of the overlap to reduce the occurences
     * of blending artifacts. If there is no overlap, the image present
     * has no blending activity, this is determined by the black
     * background and is not through an alpha layer to help reduce
     * storage needs. The algorithm below is complicated most by
     * the blending determinations, overlapping a left untransformed
     * image with a right transformed image with a black background is
     * all that remains.
     */
    /*
     * Normalize transformation against origin, the
     * transformation algorithm was in reference to the right
     * hand side of the left hand image before.
     */
    for (row = 0; row < Out->image->pam.height; row++) {
        unsigned column, firstRight = 0;
        tuple * Row = (*(Out->Row))(Out,row);
        for (column = 0; column < Out->image->pam.width; ++column) {
            unsigned plane;
            int      right_row = -1, right_column;
            right_column = (*(Stitch->XRight))(Stitch,column,
              (y < 0) ? (row + y) : row);
            if ((0 <= right_column)
             && (right_column < Right->pam.width)) {
                right_row = (*(Stitch->YRight))(Stitch,column,
                  (y < 0) ? (row + y) : row);
                if ((right_row < 0)
                 || (Right->pam.height <= right_row)) {
                    right_column = -1;
                    right_row    = -1;
                }
            } else {
                right_column = -1;
            }
            for (plane = 0; plane < Out->image->pam.depth; ++plane) {
                sample leftPixel, rightPixel;
                /* Left `mix' is easy to find */
                leftPixel = (column < Left->pam.width)
                  ? (y < 0)
                    ? ((row < -y) || (row >= (Left->pam.height - y)))
                      ? 0
                      : Left->tuple[row + y][column][plane]
                    : (row < Left->pam.height)
                      ? Left->tuple[row][column][plane]
                      : 0
                  : 0;
                rightPixel = 0;
                if (right_column >= 0) {
                    rightPixel = Right->tuple[right_row][right_column][plane];
                    if ((rightPixel > 0)
                     && (firstRight == 0)) {
                        firstRight = column;
                    }
                }
                if (leftPixel == 0) {
                    leftPixel = rightPixel;
                } else if ((firstRight <= column)
                 && (column < Left->pam.width)
                 && (rightPixel > 0)) {
                    /* blend 7/8 over half of stitch */
                    int w = Left->pam.width - firstRight;
                    if (column < (firstRight + w/2)) {
                        w = (w * 4) / 7;
                        leftPixel = (sample)(
                            ((leftPixel
                              * (unsigned long)(
                            firstRight + w
                                - column))
                            + (rightPixel
                              * (unsigned long)(column
                                - firstRight)))
                          / (unsigned long)w);
                    } else {
                        w *= 4;
                        leftPixel = (sample)(
                            ((leftPixel
                              * (unsigned long)(
                            Left->pam.width
                                - column))
                            + (rightPixel
                              * (unsigned long)(column
                            - Left->pam.width + w)))
                          / (unsigned long)w);
                    }
                }
                    
                Row[column][plane] = leftPixel;
            }
        }
        (*(Out->FlushRow))(Out,row);
    }
    (*(Out->FlushImage))(Out);
    (*(Out->DeAlloc))(Out);
    free_image (Left);
    free_image (Right);

    return 0;
} /* pnmstitch() - end */

/* Output Methods */

/* LineAtATime output method */

static void
LineAtATimeDeAlloc(Output * me)
{
    if (me->image != (Image *)NULL) {
        /* Trick the proper freeing of resouces on the Output Image */
        me->image->pam.height = 1;
        /* Free up resources */
        free_image (me->image);
        me->image = (Image *)NULL;
    }
    if (me->extra != (void *)NULL) {
        free (me->extra);
        me->extra = (void *)NULL;
    }
} /* LineAtATimeDeAlloc() - end */

static bool
LineAtATimeAlloc(Output * me, char * file, unsigned width, unsigned height, 
                 struct pam * prototype)
{
    /* Output the combined images */
    me->extra = (void *)NULL;
    me->image = allocate_image();
    if (me->image == (Image *)NULL) {
        return false;
    }
    me->image->pam = *prototype;
    me->image->pam.width = width;
    me->image->pam.height = height;
    /* Give the output a name */
    me->image->name = strdup(file);
    /* Initialize output arrays */
    if (writeinit(me->image) == false) {
        LineAtATimeDeAlloc(me);
        return false;
    }
    me->image->pam.height = 1;
    me->image->tuple = pnm_allocpamarray(&me->image->pam);
    me->image->pam.height = height;
    return true;
} /* LineAtATimeAlloc() - end */

static tuple *
LineAtATimeRow(Output * me, unsigned row)
{
    UNREFERENCED_PARAMETER(row);
    return me->image->tuple[0];
} /* LineAtATimeRow() - end */

static void
LineAtATimeFlushRow(Output * me, unsigned row)
{
    UNREFERENCED_PARAMETER(row);
    if (me->image != (Image *)NULL) {
        pnm_writepamrow(&me->image->pam, me->image->tuple[0]);
    }
} /* LineAtATimeFlushRow() - end */

static void
LineAtATimeFlushImage(Output * me)
{
    UNREFERENCED_PARAMETER(me);
} /* LineAtATimeFlushImage() - end */

/* Horizontal Crop output method */

#define HorizontalCropDeAlloc LineAtATimeDeAlloc

typedef struct {
    int state;
    int lostInSpace;
} HorizontalCropExtra;

static bool
HorizontalCropAlloc(Output * me, char * file, unsigned width, unsigned height,
                    struct pam * prototype)
{
    unsigned long pos;

    if (LineAtATimeAlloc(me, file, width, height, prototype) == false) {
        return false;
    }
    me->extra = (void *)malloc(sizeof(HorizontalCropExtra));
    if (me->extra == (void *)NULL) {
        HorizontalCropDeAlloc(me);
        return false;
    }
    memset (me->extra, 0, sizeof(HorizontalCropExtra));
    pos = ftell(me->image->pam.file);
    if ((fseek(me->image->pam.file, 1L, SEEK_SET) != 0)
     || (ftell(me->image->pam.file) != 1L)) {
        fprintf (stderr, "%s needs to output to a seekable entity\n",
          me->Name);
    }
    (void)fseek(me->image->pam.file, pos, SEEK_SET);
    return true;
} /* HorizontalCropAlloc() - end */

#define HorizontalCropRow LineAtATimeRow

static void
HorizontalCropFlushRow(Output * me, unsigned row)
{
    unsigned column;
    unsigned threshold;
    UNREFERENCED_PARAMETER(row);

    if (me->image == (Image *)NULL) {
        return;
    }
    if (((HorizontalCropExtra *)(me->extra))->state == 2) {
        ((HorizontalCropExtra *)(me->extra))->lostInSpace++;
        return;
    }
    /* Any pitch black pixels? */
    threshold = 2;
    for (column = 0; column < me->image->pam.width; ++column) {
        unsigned plane;
        for (plane = 0; 
             me->image->tuple[0][column][plane] == (sample)0; 
             ++plane) {
            if ((plane >= me->image->pam.depth)
             && (--threshold == 0)) {
                if (((HorizontalCropExtra *)(me->extra))->state == 1) {
                    ((HorizontalCropExtra *)(me->extra))->state = 2;
                }
                ((HorizontalCropExtra *)(me->extra))->lostInSpace++;
                return;
            }
        }
    }
    ((HorizontalCropExtra *)(me->extra))->state = 1;
    pnm_writepamrow(&me->image->pam, me->image->tuple[0]);
} /* HorizontalCropFlushRow() - end */

static void
HorizontalCropFlushImage(Output * me)
{
    me->image->pam.height -= ((HorizontalCropExtra *)(me->extra))->lostInSpace;
    if (verbose) {
        fprintf (stderr, "%s has set image size to %d x %d\n",
          me->Name, me->image->pam.width, me->image->pam.height);
    }
    if (fseek(me->image->pam.file, 0L, SEEK_SET) == 0) {
        pnm_writepaminit(&me->image->pam);
    } else {
        fprintf (stderr,
          "%s failed to seek to beginning to rewrite the header\n",
          me->Name);
    }
} /* HorizontalCropFlushImage() - end */

/* Output Method Table */

Output OutputMethods[] = {
    { "LineAtATime", LineAtATimeAlloc, LineAtATimeDeAlloc, LineAtATimeRow,
      LineAtATimeFlushRow, LineAtATimeFlushImage },
    { "HorizontalCrop", HorizontalCropAlloc, HorizontalCropDeAlloc,
      HorizontalCropRow, HorizontalCropFlushRow, HorizontalCropFlushImage },
    { NULL }
};

/* Stitcher Methods */

/* Linear Stitcher Methods */

static void
LinearDeAlloc(Stitcher * me)
{
    if (me->parms != (float *)NULL) {
        free (me->parms);
        me->parms = (float *)NULL;
    }
} /* LinearDeAlloc() - end */

static bool
LinearAlloc(Stitcher * me)
{
    me->parms = (float *)malloc (sizeof(float)*8);
    if (me->parms == (float *)NULL) {
        return false;
    }
    /* Constraints unset */
    me->x = INTPTR_MAX;
    me->y = INTPTR_MAX;
    me->width = INTPTR_MAX;
    me->height = INTPTR_MAX;
    return true;
} /* LinearAlloc() - end */

static void
LinearConstrain(Stitcher * me, int x, int y, int width, int height)
{
    me->x = x;
    me->y = y;
    me->width = width;
    me->height = height;
} /* LinearConstrain() - end */

/*
 *  First pass is to find an approximate match. To do so, we take a
 *  width sliver of the left hand side of the right image and compare
 *  the sample to the left hand image. Accuracy is honored over speed.
 *  The image overlap is expected between 7/16 to 1/16 in the horizontal
 *  position, and a minumum of 5/8 in the vertical dimension.
 *
 *  Blind alleys:
 *      - reduced resolution can match in totally wrong regions,
 *        as such it can not be used to improve the speed by
 *        getting close, then fine tuning at full resolution.
 *      - vector (color) average of sample matched to running
 *        vector average on left image in an attempt to improve
 *        positional accuracy of a reduced resolution image
 *        produced even more artifacts.
 *      - A complete boxed sliver did not find a minima, as for
 *        too large or too small of a square sample. heuristics
 *        show that it works between 1/128 to 1/16 of the total
 *        image dimension. Smaller, of course, improves speed.
 *
 *  Transformation parameters
 *      x=x.+a
 *      y=y'+b
 *  Where x,y represents the original point, and x.,y.
 *  represents the transformed point. Thus:
 *
 * Transformed image:
 * ((x'+x")/2,(y'+y"-H)/2)               ((x'+x"+2W)/2,(y'+y"-H)/2)
 * ((x'+x")/2,(y'+y"+H)/2)               ((x'+x"+2W)/2,(y'+y"+H)/2
 *
 * Corresponding to Original (dot) image:
 * (0,0)                 (Right->pam.width,0)
 * (0,Right->pam.height) (Right->pam.width,Right->pam.height)
 *
 *  Our matching data points are centered on x.=width/2, and
 * scan for transformation results with a variety of y. values:
 *  x=a*width/2+by.+c*width*y./2+d
 *  y=e*width/2+fy.+g*width*y./2+h
 * we set:
 *  A=b+c*width/2
 *  B=a*width/2+d
 *  C=f+g*width/2
 *  D=e*width/2+h
 * thus simplifying to:
 *  x=Ay.+B
 *  y=Cy.+D
 * adding in a weighting factor of w, the error equation is:
 *   2                       2           2
 *  E(A,B,C,D)=w * ((Ay.+B-x) + (Cy.+D-y))
 * thus
 *    2       
 *  dE(A)=2wy.(Ay.+B-x) => 0=A{wy.y. + B{wy. - {wy.x
 *    2
 *      dE(B)=2w(Ay.+B-x)   => 0=A{wy.   + B{w   - {wx
 *      A=({wy.x{w-{wx{wy.)/({wy.y.{w-{wy.{wy.)
 *      B=({wx-A{wy.)/{w
 * and
 *    2
 *      dE(C)=2wy.(Cy.+D-y) => 0=C{wy.y. + D{wy. - {wy.y
 *    2
 *      dE(D)=2w(Cy.+D-y)   => 0=C{wy.   + D{w   - {wy
 *      C=({wy.y{w-{wy{wy.)/({wy.y.{w-{wy.{wy.)
 *      D=({wy-C{wy.)/{w
 * requiring us to collect:
 *   {wy.x=sumydotx
 *     {wx=sumx
 *        {wy.=sumydot
 *      {wy.y.=sumydotydot
 *          {w=sum
 *       {wy.y=sumydoty
 *         {wy=sumy
 * Once we have A, B, C and D, we can calculate the x',y' and x",y"
 * values as follows (based on geometric interpolation from the above
 * constraints):
 *  x'=AH/2+B-AHW/(2W-width)
 *  y'=CH/2+D-CHW/(2W-width)
 *  x"=AH/2+B+AHW/(2W-width)
 *  y"=CH/2+D+CHW/(2W-width)
 * These two points can be used either in the Linear or the BiLinear to
 * establish a transform.
 */

#define IMAGE_PORTION 128
#define SKIP_SLIVER   1

#define Sliver_A   0
#define Sliver_B   1
#define Sliver_C   2
#define Sliver_D   3
#define Sliver_xp  4
#define Sliver_yp  5
#define Sliver_xpp 6
#define Sliver_ypp 7

/*
 *  Find the weighted best line fit using the left hand margin of the
 * right hand image.
 */
static bool
SliverMatch(Stitcher * me, Image * Left, Image * Right, 
            unsigned image_portion, unsigned skip_sliver)
{
    /* up/down 3/10, make sure has an odd number of members */
    int const bestSize = (1 + 2 * ((image_portion * 3) / (10 * skip_sliver)));
    Best       ** best;  /* malloc'ed array of Best *  */
    float         sumydotx, sumx, sumydot, sum;
    float         sumydoty, sumy, sumydotydot;
    int           yDiff;
    unsigned      X, Xmin, Xmax, num, xmin, xmax;
    int           x, y, Y, Ymin, Ymax, in, ymin, ymax;
    unsigned long XYcounter;
    unsigned long XYinit;

    best = malloc(bestSize * sizeof(Best *));

    /* Harry Sticks Geeses */
    if (me->width == INTPTR_MAX) {
        me->width = Right->pam.width / image_portion;
    }
    if ((me->width > (Right->pam.width/2))
        || (me->width > (Left->pam.width/2))) {
        pm_error ("stitch sample too wide %d\n", me->width);
    } 
    if (me->height == INTPTR_MAX) {
        me->height = Right->pam.height / image_portion;
    }
    if ((me->height > Right->pam.height)
        || (me->height > Left->pam.height)) {
        pm_error ("stitch sample too high %d\n", me->height);
    }
    /* Calculate overlap position, if not overriden */
    if (me->x == INTPTR_MAX || me->y == INTPTR_MAX) {
        yDiff = (Right->pam.height * skip_sliver) / image_portion;
        XYinit = XYcounter = (unsigned long)-1L;

        /* Allocate */
        for (X = 0; X < bestSize; ++X) {
            best[X] = allocate_best();
        }
        /* Constrain loop */
        if (me->x == INTPTR_MAX) {
            Xmax = Left->pam.width - me->width;
            /* I can't bring myself to go half way */
            Xmin = Left->pam.width - (7 * Right->pam.width / 16);
        } else {
            Xmin = Xmax = me->x;
            ++Xmax;
        }
        if (me->y == INTPTR_MAX) {
            /* Middle 1/4 */
            Ymin = Y = Left->pam.height * 3 / 8;
            Ymax = Left->pam.height - Y - me->height;
        } else {
            Ymax = Ymin = me->y;
            ++Ymax;
        }
        /* Find the best */
        if (verbose) {
            fprintf (stderr, "Test %d<x<%d %d<y<%d\n",
                     Xmin, Xmax, Ymin, Ymax);
            if ((verbose == 1) || (verbose == 2)) {
                fprintf (stderr, "%79s|\r|", "");
                XYinit = (unsigned long)((
                    (unsigned long)(Xmax - Xmin)
                    * (unsigned long)(Ymax - Ymin))
                    * (unsigned long) bestSize 
                    / 78L);
                XYcounter = XYinit;
            }
        }
        /* A point in the middle of the right image */
        x = 0;
        y = (Right->pam.height - me->height) / 2;
        /*
         *  Exhaustively search for the best match, improvements
         * in the algorithm here, if any, would give us the best
         * bang for the buck when it comes to improving performance.
         */
        /* Loop through the constraints to find the best match */
        sumydotx=sumx=sumydot=sumydotydot=sum=sumydoty=sumy=0.0;
        xmin = UINTPTR_MAX;
        xmax = 0;
        ymin = INTPTR_MAX;
        ymax = INTPTR_MIN;
        for (in = num = 0;; ++num) {
            float w;
            int offY = yDiff * (in - bestSize/2);
            for (X = Xmin; X < Xmax; ++X) {
                int theX = X - Left->pam.width;
                for (Y = Ymin; Y < Ymax; ++Y) {
                    unsigned long total = sub_tuple(
                        Left, X, Y + offY,
                        Right, x, y + offY,
                        me->width, me->height);
                    update_best(best[in], total, theX, Y + offY);
                    if (--XYcounter == 0) {
                        XYcounter = XYinit;
                        fprintf (stderr, "*");
                    }
                }
            }
            /* slop */
            {
                float SUMx, SUMxx, SUMy, SUMyy, SUMw;
                SUMx = SUMxx = SUMy = SUMyy = SUMw = 0.0;
                for (X = 0; X < NUM_BEST; ++X) {
                    float const w2 = best[in][X].total > 0 ?
                        1.0 / (float)best[in][X].total : 1.0;
                    SUMx  += w2 * best[in][X].x;
                    SUMy  += w2 * best[in][X].y;
                    SUMxx += w2 * (best[in][X].x * best[in][X].x);
                    SUMyy += w2 * (best[in][X].y * best[in][X].y);
                    SUMw  += w2;
                }
                /* Find our weighted error */
                w = SUMw / 
                    ((SUMxx - (SUMx*SUMx)/SUMw) + (SUMyy - (SUMy*SUMy)/SUMw));
            }
            /* magnify slop */
            w *= w;
            Y = y + offY;
            sumy        += w * best[in][0].y;
            sumx        += w * best[in][0].x;
            sum         += w;
            sumydot     += w * Y;
            sumydotydot += w * Y * Y;
            sumydoty    += w * Y * best[in][0].y;
            sumydotx    += w * Y * best[in][0].x;
            /* Calculate the best fit line for these matches */
            me->parms[Sliver_C] = ((sumydotydot * sum)
                                   - (sumydot * sumydot));
            if (me->parms[Sliver_C] == 0.0) {
                me->parms[Sliver_A] = 0.0;
            } else {
                me->parms[Sliver_A] = ((sumydotx * sum)
                                       - (sumx * sumydot))
                    / me->parms[Sliver_C];
                me->parms[Sliver_C] = ((sumydoty * sum)
                                       - (sumy * sumydot))
                    / me->parms[Sliver_C];
            }
            if (sum == 0.0) {
                me->parms[Sliver_B] = me->parms[Sliver_D] = 0;
            } else {
                me->parms[Sliver_B] = (sumx
                                       - (me->parms[Sliver_A] * sumydot))
                    / sum;
                me->parms[Sliver_D] = (sumy
                                       - (me->parms[Sliver_C] * sumydot))
                    / sum;
            }
            if (verbose > 2) {
                fprintf (stderr, "%.4g*(%d,%d)@(%d,%d)\n",
                         w, best[in][0].x + Left->pam.width, best[in][0].y,
                         me->width / 2, Y);
            }
            /* Record history of limits */
            if ((best[in][0].x + Left->pam.width) < xmin) {
                xmin = best[in][0].x + Left->pam.width;
            }
            if (xmax < (best[in][0].x + Left->pam.width)) {
                xmax = best[in][0].x + Left->pam.width;
            }
            if ((best[in][0].y - offY) < ymin) {
                ymin = best[in][0].y - offY;
            }
            if (ymax < (best[in][0].y - offY)) {
                ymax = best[in][0].y - offY;
            }
            /* Lets restrict the search a bit now */
            if (num > 0) {
                if (me->x == INTPTR_MAX) {
                    int newXmin, newXmax, hold;
                    /* Use the formula to determine the bounds */
                    newXmin = (int)(me->parms[Sliver_B] + 0.5)
                        + Left->pam.width;
                    newXmax = (int)((me->parms[Sliver_A]
                                     * (float)Left->pam.height)
                                    + me->parms[Sliver_B] + 0.5)
                        + Left->pam.width;
                    if (newXmax < newXmin) {
                        hold = newXmin;
                        newXmin = newXmax;
                        newXmax = hold;
                    }
                    /* Trust little ... */
                    hold = (3 * newXmin - newXmax) / 2;
                    newXmax = (3 * newXmax - newXmin) / 2;
                    newXmin = hold;
                    /* Don't go inside history */
                    if (xmin < newXmin) {
                        newXmin = xmin + Left->pam.width;
                    }
                    if (newXmax < xmax) {
                        newXmax = xmax + Left->pam.width;
                    }
                    /* Now upgrade new minimum and maximum */
                    hold = Xmin;
                    if ((Xmin < newXmin) && (newXmin < Xmax)) {
                        hold = newXmin;
                    }
                    if ((Xmin < newXmax) && (newXmax < Xmax)) {
                        Xmax = newXmax;
                    }
                    Xmin = hold;
                }
                if (me->y == INTPTR_MAX) {
                    int newYmin, newYmax, hold;
                    float tmp;
                    /* Use the formula to determine the bounds */
                    newYmin = tmp = me->parms[Sliver_D]
                        + ((float)(Left->pam.height + 1))
                        / 2;
                    newYmax = (int)((me->parms[Sliver_C]
                                     * (float)Left->pam.height)
                                    + tmp) - Left->pam.height;
                    if (newYmax < newYmin) {
                        hold = newYmin;
                        newYmin = newYmax;
                        newYmax = hold;
                    }
                    /* Trust little ... */
                    hold = (3 * newYmin - newYmax) / 2;
                    newYmax = (3 * newYmax - newYmin) / 2;
                    newYmin = hold;
                    /* Don't go inside history */
                    if (ymin < newYmin) {
                        newYmin = ymin;
                    }
                    if (newYmax < ymax) {
                        newYmax = ymax;
                    }
                    /* Now upgrade new minimum and maximum */
                    hold = Ymin;
                    if ((Ymin < newYmin) && (newYmin < Ymax)) {
                        hold = newYmin;
                    }
                    if ((Ymin < newYmax) && (newYmax < Ymax)) {
                        Ymax = newYmax;
                    }
                    Ymin = hold;
                }
                if ((verbose == 1) || (verbose == 2)) {
                    XYinit = (unsigned long)(
                        (
                            (unsigned long)(Xmax - Xmin)
                            * (unsigned long)(Ymax - Ymin))
                        * (unsigned long) bestSize
                        / 78L);
                    if (XYcounter > XYinit) {
                        XYcounter = XYinit;
                    }
                }
            }
            if (in > bestSize/2) {
                in = bestSize - in;
            } else if (in < bestSize/2) {
                in = (bestSize-1) - in;
            } else {
                break;
            }
        }
        if ((verbose == 1) || (verbose == 2)) {
            fprintf (stderr, "\n");
        }
        if (verbose > 2) {
            fprintf (stderr, "Up  ");
            pr_best(best[bestSize-1]);
            fprintf (stderr, "\nMid ");
            pr_best(best[bestSize/2]);
            fprintf (stderr, "\nDown");
            pr_best(best[0]);
            fprintf (stderr, "\n");
        }

        if (verbose) {
            if (verbose > 1) {
                fprintf (stderr,
                         "[y=%.4g [x=%.4g [=%.4g [y.=%.4g "
                         "[y.y.=%.4g [y.y=%.4g [y.x=%.4g\n",
                         sumy, sumx, sum, sumydot, sumydotydot,
                         sumydoty, sumydotx);
            }
            fprintf (stderr, "x=%.4gY%+.4g\ny=%.4gY%+.4g\n",
                     me->parms[Sliver_A], me->parms[Sliver_B],
                     me->parms[Sliver_C], me->parms[Sliver_D]);
        }
        /*
     *  Free up resources
     */
        for (X = 0; X < bestSize; ++X) {
            free (best[X]);
        }
        free(best);

        /* Calculate x',y' and x",y" from best fit line formula */
        sum = (float)(Right->pam.width * Right->pam.height)
            / (float)(2 * Right->pam.width - me->width);
        me->parms[Sliver_xpp] = me->parms[Sliver_A] * sum;
        me->parms[Sliver_ypp] = me->parms[Sliver_C] * sum;
        sumx = me->parms[Sliver_A] * (Right->pam.height / 2)
            + me->parms[Sliver_B];
        sumx += Left->pam.width;
        sumy = me->parms[Sliver_C] * (Right->pam.height / 2)
            + me->parms[Sliver_D];
        me->parms[Sliver_xp] = sumx - me->parms[Sliver_xpp];
        me->parms[Sliver_yp] = sumy - me->parms[Sliver_ypp];
        me->parms[Sliver_xpp] += sumx;
        me->parms[Sliver_ypp] += sumy;
        if (verbose > 1) {
            fprintf (stderr, "x',y'=%.4g,%.4g x\",y\"=%.4g,%.4g\n",
                     me->parms[Sliver_xp], me->parms[Sliver_yp],
                     me->parms[Sliver_xpp], me->parms[Sliver_ypp]);
        }
    }
    return true;
} /* SliverMatch() - end */

#define Linear_a   0
#define Linear_b   1
#define Linear_c   2
#define Linear_d   3
#define Linear_e   4
#define Linear_f   5
#define Linear_g   6
#define Linear_h   7

static bool
LinearMatch(Stitcher * me, Image * Left, Image * Right)
{
    if (SliverMatch(me, Left, Right, IMAGE_PORTION, SKIP_SLIVER * 8) == 
        false) {
        return false;
    }

    me->x = - (me->parms[Sliver_xp] + me->parms[Sliver_xpp] + 1) / 2;
    me->y = - (me->parms[Sliver_yp] + me->parms[Sliver_ypp] + 
               (1 - Left->pam.height)) / 2;

    return true;
} /* LinearMatch() - end */

/*
 *  Transformation parameters
 *      x`=x
 *      y`=y
 */
static int
LinearXLeft(Stitcher * me, int x, int y)
{
    UNREFERENCED_PARAMETER(y);
    return x;
} /* LinearXLeft() - end */

static int
LinearYLeft(Stitcher * me, int x, int y)
{
    UNREFERENCED_PARAMETER(x);
    return y;
} /* LinearYLeft() - end */

/*
 *  Transformation parameters
 *      x`=x+me->x
 *      y`=y+me->y
 */
static int
LinearXRight(Stitcher * me, int x, int y)
{
    UNREFERENCED_PARAMETER(y);
    return x + me->x;
} /* LinearXRight() - end */

static int
LinearYRight(Stitcher * me, int x, int y)
{
    UNREFERENCED_PARAMETER(x);
    return y + me->y;
} /* LinearYLeft() - end */

static void
LinearOutput(Stitcher * me, FILE * fp)
{
    fprintf (fp, "x'=x%+d\ny'=y%+d\n", me->x, me->y);
} /* LinearOutput() - end */

/* BiLinear Stitcher Methods */

#define BiLinearDeAlloc LinearDeAlloc

static bool
BiLinearAlloc(Stitcher * me)
{
    unsigned num;
    if (LinearAlloc(me) == false) {
        return false;
    }
    for (num = 8; num > 0; --num) {
        me->parms[num] = 0.0;
    }
    /* Unity Transform matrix */
    me->parms[0] = me->parms[5] = 1.0;
    return true;
} /* BiLinearAlloc() - end */

static void
BiLinearConstrain(Stitcher * me, int x, int y, int width, int height)
{
    LinearConstrain(me, x, y, width, height);
    if (x != INTPTR_MAX) {
        me->parms[3] -= x;
    }
    if (y != INTPTR_MAX) {
        me->parms[7] -= y;
    }
} /* BiLinearConstrain() - end */

/*
 *  First pass is to find an approximate match. To do so, we take a
 *  width sliver of the left hand side of the right image and compare
 *  the sample to the left hand image. Accuracy is honored over speed.
 *  The image overlap is expected between 7/16 to 1/16 in the horizontal
 *  position, and a minumum of 5/8 in the vertical dimension.
 *
 *  Blind alleys:
 *      - Tried a simpler constraint for right side to be `back'
 *        to image, twisted too much sometimes:
 *         . . .
 *         W=aW+bH+cWH+d
 *         H=eW+fH+gWH+h
 *         W=aW+d
 *         0=eW+h
 *        Solve the equations resulted in:
 *         a = W/(W-x") - cy"
 *         b = -Wc
 *         c = W/((x"-W)(y'-y"))
 *         d = (1-a)W
 *         e = y'(y"x'+W-Hx'-x"y")/(x'y"x"-Wx'y"-Wy'x"-WWy'+x'y'x"-Wx'y'-W)
 *         f = 1 - Wg
 *         g = (e + (H-y")/(x"-W))/y"
 *         h = -We
 *        Results left here for historical reasons.
 *
 *  Transformation parameters
 *      x=ax.+by.+cx.y.+d
 *      y=ex.+fy.+gx.y.+h
 *  Where x,y represents the original point, and x.,y.
 *  represents the transformed point. Thus:
 *
 * Transformed image:
 * (x',y')               (x'",y'")
 * (x",y")               (x"",y"")
 *
 * Corresponding to Original (dot) image:
 * (0,0)                 (Right->pam.width,0)
 * (0,Right->pam.height) (Right->pam.width,Right->pam.height)
 *
 * Define:
 *      H=Right->pam.height
 *      W=Right->pam.width
 * Given that I want a flat presentation that both reduces the distortion
 * necessary on an image, reduces the cropping losses, and flattens out the
 * spherical or orbit distortions; it was chosen to constrain the right side
 * in the middle horizontal, and pivot the left side in that middle (hopefully
 * minimally) and to allow the image only vertical and horizontal location
 * placement. Rotating the entire image could increase cropping losses
 * especially if the focus was not down the center of the image on a
 * graduated field causing the distortion to accumulate in subsequent
 * images. Trapezoidal would cause the distortion to accumulate in subsequent
 * images as well, resetting to `square' gradually towards the right would
 * allow the next image to restart a match placing the distortions mainly
 * in the stitching zone where averaging and the slight expectation of
 * artifacts would minimize the effects. These constraints can be explained
 * mathematically as the following:
 *  x'" + x"" - 2W = x' + x"
 *       y'" + y"" = y' + y"
 *                 x'" = x""
 *         y'" + H = y""
 * resulting in the right side of the image being completely explained by the
 * placement of the left hand side:
 *  x'"=(x'+x"+2W)/2
 *  y'"=(y'+y"-H)/2
 *  x""=(x'+x"+2W)/2
 *  y""=(y'+y"+H)/2
 *
 * Describing the `X' polygon using geometry and ratios:
 *  X=A(y-(y'+y")/2)(x-(x'+x"+2W)/2) + (x - (x'+x")/2))
 *    A=2(x"-x')/((y'-y")(x'-x"-2W))
 *  a=2(y'(x'-x")-W(y'-y"))/((y'-y")(x'-x"-2W))
 *  b=(x'-x")(x'+x"+2W)/((y'-y")(x'-x"-2W))
 *  c=2(x"-x')/((y'-y")(x'-x"-2W))
 *  d=(2W(x"y'-x'y")+y'(x"x"-x'x'))/((y'-y")(x'-x"-2W))
 *
 * Describing the `Y' polygon using geometry and ratios (note use of X rather
 * than x, this has the effect of linearalizing the polygon).
 *  Y=((y'-y"+H)/W(y'-y"))(y-(y'+y")/2)(X-HW/(y'-y"+H)) + H/2
 *  e=(y'+y")(y'-y"+H)/2W(y"-y')
 *  f=H/(y"-y')
 *  g=(y'-y"+H)/W(y'-y")
 *  h=Hy'/(y'-y")
 *
 * FYI: Reverse transform using the same formula style is:
 *  a=(x"-x'+2W)/2W
 *  b=(x"-x')/H
 *  c=(x'-x")/WH
 *  d=x'
 *  e=(y"-y'-H)/2W
 *  f=(y"-y')/H
 *  g=(y'-y"+H)/WH
 *  h=y'
 */

#define BiLinear_a   0
#define BiLinear_b   1
#define BiLinear_c   2
#define BiLinear_d   3
#define BiLinear_e   4
#define BiLinear_f   5
#define BiLinear_g   6
#define BiLinear_h   7

static bool
BiLinearMatch(Stitcher * me, Image * Left, Image * Right)
{
    /* float A, B, C, D; */
    float xp, yp, xpp, ypp;

    if (SliverMatch(me, Left, Right, IMAGE_PORTION, SKIP_SLIVER) == false) {
        return false;
    }
    /*
    A   = me->parms[Sliver_A];
    B   = me->parms[Sliver_B];
    C   = me->parms[Sliver_C];
    D   = me->parms[Sliver_D];
     */
    xp  = me->parms[Sliver_xp];
    yp  = me->parms[Sliver_yp];
    xpp = me->parms[Sliver_xpp];
    ypp = me->parms[Sliver_ypp];

    /*
     *  Calculate any necessary transformations on the
     * right image to improve the stitching match. We have Done a
     * weighted best fit line on the points we have collected
     * thus far, now translate this to the constrained
     * transformation equations.
     */
    /* a = y"-y' */
    me->parms[BiLinear_a] = ypp-yp;
    /* c = x'-x" */
    me->parms[BiLinear_c] = xp-xpp;
    /* d = (y"-y')(x"-x'+2W) = (y'-y")(x'-x"-2W) */
    me->parms[BiLinear_d] = me->parms[BiLinear_a]
                            * ((float)
                               (2*Right->pam.width)-me->parms[BiLinear_c]);
    /* a = 2(y'(x'-x")+W(y"-y'))/((y'-y")(x'-x"-2W)) */
    me->parms[BiLinear_a] = 2*(yp*me->parms[BiLinear_c]
                          + me->parms[BiLinear_a]*(float)(Right->pam.width))
                          / me->parms[BiLinear_d];
    /* b = (x'-x")(x'+x"+2W)/((y'-y")(x'-x"-2W)) */
    me->parms[BiLinear_b] = me->parms[BiLinear_c]
                          * (xp+xpp+(float)(2*Right->pam.width))
                          / me->parms[BiLinear_d];
    /* c = -2(x'-x")/((y'-y")(x'-x"-2W)) */
    me->parms[BiLinear_c]*= -2/me->parms[BiLinear_d];
    /* d = (2W(x"y'-x'y")+y'(x"x"-x'x'))/((y'-y")(x'-x"-2W)) */
    me->parms[BiLinear_d] = ((xpp*yp-xp*ypp)*(float)(2*Right->pam.width)
                          + yp*(xpp*xpp-xp*xp))
                          / me->parms[BiLinear_d];

    /* f = y"-y' */
    me->parms[BiLinear_f] = ypp-yp;
    /* g = (y"-y'-H)/W(y"-y') */
    me->parms[BiLinear_g] = (me->parms[BiLinear_f]-(float)Right->pam.height)
                          / me->parms[BiLinear_f]
                          / (float)Right->pam.width;
    /* e = (y'+y")(y'-y"+H)/2W(y"-y') = -g(y'+y")/2 */
    me->parms[BiLinear_e] = (yp+ypp)*me->parms[BiLinear_g]/-2;
    /* f=H/(y"-y') */
    me->parms[BiLinear_f] = ((float)Right->pam.height)
                          / me->parms[BiLinear_f];
    /* h = Hy'/(y'-y") = -fy' */
    me->parms[BiLinear_h] = -yp*me->parms[BiLinear_f];

    return true;
} /* BiLinearMatch() - end */

/*
 *  Transformation parameters
 *      x`=x
 *      y`=y
 */
#define BiLinearXLeft LinearXLeft
#define BiLinearYLeft LinearYLeft

/*
 *  Transformation parameters
 *      x`=ax+by+cxy+d
 *      y`=ex`+fy+gx`y+h
 */
static int
BiLinearXRight(Stitcher * me, int x, int y)
{
    return (me->parms[BiLinear_a] * x) + (me->parms[BiLinear_b] * y)
         + (me->parms[BiLinear_c] * (x * y)) + me->parms[BiLinear_d] + 0.5;
} /* BiLinearXRight() - end */

static int
BiLinearYRight(Stitcher * me, int x, int y)
{
    /* A little trick I learned from a biker */
    x = BiLinearXRight(me, x, y);
    return (me->parms[BiLinear_e] * x) + (me->parms[BiLinear_f] * y)
         + (me->parms[BiLinear_g] * (x * y)) + me->parms[BiLinear_h] + 0.5;
} /* BiLinearYLeft() - end */

static void
BiLinearOutput(Stitcher * me, FILE * fp)
{
    fprintf (fp,
      "x'=%.6gx%+.6gy%+.6gxy%+.6g\ny'=%.6gx'%+.6gy%+.6gx'y%+.6g\n",
      me->parms[0], me->parms[1], me->parms[2], me->parms[3],
      me->parms[4], me->parms[5], me->parms[6], me->parms[7]);
} /* BiLinearOutput() - end */

/* Stitcher Method Table */

Stitcher StitcherMethods[] = {
    { "BiLinearSliver", BiLinearAlloc, BiLinearDeAlloc, BiLinearConstrain,
      BiLinearMatch, BiLinearXLeft, BiLinearYLeft, BiLinearXRight,
      BiLinearYRight, BiLinearOutput },
    { "LinearSliver", LinearAlloc, LinearDeAlloc, LinearConstrain,
      LinearMatch, LinearXLeft, LinearYLeft, LinearXRight,
      LinearYRight, LinearOutput },
    { NULL }
};

