/*
* 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 software, even if
* advised of the possibility of such damage.
*
* Any restrictions or encumberances added to this source code or derivitives,
* is prohibited.
*
* Name: pnmstitch.c
* Description: Automated panoramic stitcher.
* Many digital Cameras have a panorama mode where they hold on to
* the right hand side of an image, shifted to the left hand side of their
* view screen for subsequent pictures, facilitating the manual stitching
* up of a panoramic shot. These same cameras are shipped with software
* to manually or automatically stitch images together into the composite
* image. However, these programs are dedicated for a specific OS.
* In the pnmstitch program, it analyzes the match between the images,
* generates a transform, processes the transform on the images, and
* blends the overlapping regions. In addition, there is an output filter
* to process automatic cropping of the resultant image.
* The stitching software here works by constraining the right
* hand side of the right hand image as `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 loose term) of the left hand
* image.
* Author: Mark Salyzyn <mark@bohica.net> June 2002
* Version: 0.0.4
*
* Modifications: 0.0.4 July 31 2002 Mark Salyzyn <mark@bohica.net>
* & Bryan Henderson <bryanh@giraffe-data.com>
* - FreeBSD port.
* - merge changes to incorporate into netpbm tree.
* Modifications: 0.0.3 July 27 2002 Mark Salyzyn <mark@bohica.net>
* & "George M. Sipe" <geo@sipe.org>
* - Deal with subtle differences between BSD and GNU getopt
* facilitating the Linux port.
* Modifications: 0.0.2 July 25 2002 Mark Salyzyn <mark@bohica.net>
* - RotateSliver needs to use higher resolution match.
* - RotateCrop code interfered with StraightThrough code
* resulting in an incorrect pnm image output.
* Modifications: 0.0.1 July 18 2002 Mark Salyzyn <mark@bohica.net>
* - Added BiLinearSliver, RotateSliver and HorizontalCrop
*
* ToDo:
* - Split this into multiple files ... nah, keep it in one to
* keep it from polluting the netpbm tree.
* - Add and refine the videorbits algorithm. One piece of public
* domain software that is pnm aware is called videorbits. It
* needs considerably more overlap than the Digital Cameras set
* up (of course, anyone can to a panorama with more overlap
* with our without the feature) as it uses what is called video
* flow to generate the match. It is designed more for a series
* of images from a video camera. videorbits has three programs,
* one generates the transform, the next processes the
* transform, and the final blends the images together much as
* this one piece program does.
* - speedups in matching algorithm
* - refinement in accuracy of matching algorithm
* - Add RotateCrop filter algorithm (in-memory copy of image,
* detect least loss horizontal crop on a rotated image).
* - pnmstitch should be generalized to handle transformation
* occurring on the left image, currently it blends assuming
* that there is no transformation effects on the left image.
* - user selectable blending algorithms?
*/
#define _DEFAULT_SOURCE 1 /* New name for SVID & BSD source defines */
#define _BSD_SOURCE 1 /* Make sure strdup() is in string.h */
#define _XOPEN_SOURCE 500 /* Make sure strdup() is in string.h */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include "pm_c_util.h"
#include "shhopt.h"
#include "nstring.h"
#include "mallocvar.h"
#include "pam.h"
/*
* Structures
*/
/*
* Image structure
*/
typedef struct {
const char * name; /* File Name */
struct pam pam; /* netpbm image description */
tuple ** tuple; /* in-memory copy of image */
} Image;
/*
* Output class
* The following methods and data allocations are used for output filter.
*/
typedef struct output {
/* name */
const char * Name;
/* methods */
bool (* Alloc)(struct output * me, const char * file,
unsigned int width, unsigned int 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
* The following methods and data allocations are used for operations
* surrounding stitching of an image.
*/
typedef struct stitcher {
/* name */
const char * Name;
/* methods */
bool (* Alloc)(struct stitcher *me);
void (* DeAlloc)(struct stitcher *me);
void (* Constrain)(struct stitcher *me, int x, int y,
int width, int height);
/* Set transformation parameter constraints. This affects the
function of a future 'Match' method execution.
*/
bool (* Match)(struct stitcher *me, Image * Left, Image * Right);
/* Determine the transformation parameters for the stitching.
I.e. determine the parameters that affect future invocations
of the transformation methods below. You must execute a
'Match' before executing any of the transformation methods.
*/
/*-----------------------------------------------------------------------
The transformation methods answer the question, "Which pixel in the left
image and which pixel in the right image contribute to the pixel at
Column X, Column Y of the output?
If there is no pixel in the left image that contributes to the output
pixel in question, the methods return column or row numbers outside
the bounds of the left image (possibly negative). Likewise for the
right image.
*/
float (* XLeft)(struct stitcher *me, int x, int y);
/* column number of the pixel from the left image */
float (* YLeft)(struct stitcher *me, int x, int y);
/* row number of the pixel from the left image */
float (* XRight)(struct stitcher *me, int x, int y);
/* column number of the pixel from the right image */
float (* YRight)(struct stitcher *me, int x, int y);
/* row number of the pixel from the left image */
/*----------------------------------------------------------------------*/
/* Output methods */
void (* Output)(struct stitcher *me, FILE * fp);
/* private data */
int x, y, width, height;
/* For a Linear Sliver stitcher, 'x' and 'y' are simply the offset you
add to an output location to get the location in the right image of the
pixel that corresponds to that output pixel.
*/
float * parms;
} Stitcher;
extern Stitcher StitcherMethods[];
/*
* Prototypes
*/
static int pnmstitch(const char * const left,
const char * const right,
const char * const out,
int const x,
int const y,
int const width,
int const height,
const char * const stitcher,
const char * const filter);
struct cmdlineInfo {
/*
* All the information the user supplied in the command line,
* in a form easy for the program to use.
*/
const char * leftFilespec; /* '-' if stdin */
const char * rightFilespec; /* '-' if stdin */
const char * outputFilespec; /* '-' if stdout */
const char * stitcher;
const 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 pm_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 */
pm_optParseOptions3( &argc, argv, opt, sizeof(opt), 0);
/* Uses and sets argc, argv, and some of *cmdlineP and others. */
if (!widthSpec) {
cmdlineP->width = INT_MAX;
}
if (!heightSpec) {
cmdlineP->height = INT_MAX;
}
if (!xrightposSpec) {
cmdlineP->xrightpos = INT_MAX;
}
if (!yrightposSpec) {
cmdlineP->yrightpos = INT_MAX;
}
if (!stitcherSpec) {
cmdlineP->stitcher = "BiLinearSliver";
}
if (!filterSpec) {
cmdlineP->filter = "StraightThrough";
}
if (argc-1 > 3) {
pm_error("Program takes at most three arguments: left, right, and "
"output file specifications. You specified %d", argc-1);
/* NOTREACHED */
} 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.");
/* NOTREACHED */
} else if (argc-1 == 3) {
cmdlineP->outputFilespec = argv[3];
} else if (outputSpec) {
cmdlineP->outputFilespec = outputOpt;
} else {
cmdlineP->outputFilespec = "-";
}
}
} /* parseCommandLine() - end */
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;
MALLOCVAR(retVal);
if (retVal != NULL)
memset (retVal, 0, (unsigned)sizeof(Image));
return retVal;
}
/*
* free an image structure.
*/
static void
free_image(Image * image)
{
if (image->name) {
pm_strfree(image->name);
image->name = NULL;
}
if (image->tuple) {
pnm_freepamarray(image->tuple, &image->pam);
image->tuple = NULL;
}
if (image->pam.file) {
fclose (image->pam.file);
image->pam.file = NULL;
}
free (image);
}
static void
openWithPossibleExtension(const char * const baseName,
FILE ** const ifPP,
const char ** const filenameP) {
/* list of possible extensions for input file */
const char * const extlist[] = {
"", ".pnm", ".pam", ".pgm", ".pbm", ".ppm"
};
FILE * ifP;
unsigned int extIndex;
ifP = NULL; /* initial value -- no file opened yet */
for (extIndex = 0; extIndex < ARRAY_SIZE(extlist) && !ifP; ++extIndex) {
const char * trialName;
pm_asprintf(&trialName, "%s%s", baseName, extlist[extIndex]);
ifP = fopen(trialName, "rb");
if (ifP)
*filenameP = trialName;
else
pm_strfree(trialName);
}
if (!ifP)
pm_error ("Failed to open input file named '%s' "
"or '%s' with one of various common extensions.",
baseName, baseName);
*ifPP = ifP;
}
/*
* Create an image object for the PNM/PAM image in the file name 'name'
* This includes reading the entire image.
*/
static Image *
readinit(const char * const name)
{
Image * const image = allocate_image();
Image * retVal;
if (image == NULL)
retVal = NULL;
else {
FILE * ifP;
if (streq(name, "-")) {
ifP = stdin;
image->name = strdup("<stdin>");
} else {
openWithPossibleExtension(name, &ifP, &image->name);
}
image->tuple = pnm_readpam(ifP, &(image->pam),
PAM_STRUCT_SIZE(tuple_type));
fclose (ifP);
image->pam.file = NULL;
if (image->tuple == NULL) {
free_image(image);
retVal = NULL;
} else
retVal = image;
}
return retVal;
} /* readinit() - end */
/*
* Prepare an image to be output.
* Too bad we can't help them and add a .pnm on the filename,
* since this would be `bad'. on readinit we can check if the .pnm
* needs to be added ...
*/
static bool
writeinit(Image * image)
{
if (streq(image->name, "-")) {
image->pam.file = stdout;
pm_strfree(image->name);
image->name = strdup("<stdout>");
} else {
image->pam.file = pm_openw(image->name);
}
return TRUE;
} /* writeinit() - end */
/*
* Compare a subimage to an image.
* The most time consuming actions surround this subroutine.
* Return the magnitude of the difference between the specified region
* in image 'left' and the specified region in image 'right'.
* 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.
*
* The region in the left image is the rectangle with top left corner at
* Column 'lx', Row 'ly', with dimensions 'width' columns by 'height'
* rows. The region in the right image is a rectangle the same dimensions
* with upper left corner at Column 'rx', Row 'ry'.
*
* Caller must ensure that the regions indicated are entirely within the
* their respective images.
*/
static unsigned long
regionDifference(Image * left,
int lx,
int ly,
Image * right,
int rx,
int ry,
int width,
int height)
{
unsigned long total;
unsigned row;
total = 0; /* initial value */
for (row = 0; row < height; ++row) {
unsigned column;
for (column = 0; column < width; ++column) {
unsigned plane;
for (plane = 0; plane < left->pam.depth; ++plane) {
sample const leftSample = left->tuple[row][column][plane];
sample const rightSample = right->tuple[row][column][plane];
total += SQR(leftSample - rightSample);
}
}
}
return total;
}
/*
* 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
/*
* Allocate the Best structure.
*/
static Best *
allocate_best(void)
{
Best * retVal;
MALLOCARRAY(retVal, NUM_BEST);
if (retVal != NULL) {
unsigned int i;
for (i = 0; i < NUM_BEST; ++i) {
retVal[i].total = ULONG_MAX;
retVal[i].x = INT_MAX;
retVal[i].y = INT_MAX;
}
}
return retVal;
}
/*
* Free the Best structure
*/
#define free_best(best) free(best)
/*
* Placement helper for the Best structure.
*/
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;
}
/*
* Print helper for the Best structure.
*/
static void
pr_best(Best * const best)
{
int i;
if (best != (Best *)NULL)
for (i = 0; i < NUM_BEST; ++i) {
fprintf (stderr, " (%d,%d)%lu",
best[i].x, best[i].y, best[i].total);
}
} /* pr_best() - end */
static void *
findObject(const char * const name, void * start, unsigned size)
{
char ** object;
char ** best;
unsigned length;
if (name == (char *)NULL) {
return start;
}
for (length = 0, best = (char **)NULL, object = (char **)start;
*object != (char *)NULL;
object = (char **)(((char *)object) + size)) {
const char * np = name;
char * op = *object;
unsigned matched = 0;
/* case insensitive match */
while ((*np != '\0') && (*op != '\0')
&& ((*np == *op) || (tolower(*np) == tolower(*op)))) {
++np;
++op;
++matched;
}
if ((*np == '\0') && (*op == '\0')) {
break;
}
if ((matched >= length) && (*np == '\0')) {
if (matched == length) {
best = (char **)NULL;
} else {
best = object;
}
length = matched;
}
}
if (*object == (char *)NULL) {
object = best;
}
if (object == (char **)NULL) {
fprintf (stderr,
"Unknown driver \"%s\". available drivers are:\n", name);
for (object = (char **)start;
*object != (char *)NULL;
object = (char **)(((char *)object) + size)) {
fprintf (stderr, "\t%s%s\n", *object,
(object == (char **)start) ? " (default)" : "");
}
}
return object;
}
/*
* The general wrapper for both the Output and the Stitcher algorithms.
*/
/* Determine the mask corners for both Left and Right images */
static void
determineMaskCorners(Stitcher * Stitch,
Image * Left,
Image * Right,
int xp[],
int yp[])
{
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 x, y, xx, yy, count = 65536; /* max iterations */
float (*X)(Stitcher *me, int x, int y);
float (*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]) + 0.5) - 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]) + 0.5) - 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",
"es[(i%4)*3], "es[(i%4)*3],
xp[i], yp[i]);
}
}
}
} /* determineMaskCorners() - end */
static void
calculateXyWidthHeight(int xp[],
int yp[],
int * const xP,
int * const yP,
int * const widthP,
int * const heightP)
{
int x, y, width, height, i;
/* Calculate generic x,y left top corner, and the width and height */
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;
} /* calculateXyWidthHeight() - end */
static void
printPlan(int xp[], int yp[], Image * Left, Image * Right)
{
/* Calculate Left image transformed bounds */
int X, Y, W, H, 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);
} /* printPlan() - end */
static void
stitchOnePixel(Image * const Left,
Image * const Right,
struct pam const outpam,
int const row,
int const column,
int const y,
int const right_row,
int const right_column,
unsigned * const firstRightP,
tuple const outPixel) {
unsigned plane;
for (plane = 0; plane < outpam.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) && (*firstRightP == 0))
*firstRightP = column;
}
if (leftPixel == 0) {
leftPixel = rightPixel;
} else if ((*firstRightP <= column)
&& (column < Left->pam.width)
&& (rightPixel > 0)) {
/* blend 7/8 over half of stitch */
int const w = Left->pam.width - *firstRightP;
if (column < (*firstRightP + w/2)) {
int const v = (w * 4) / 7;
leftPixel = (sample)(
((leftPixel
* (unsigned long)(*firstRightP + v - column))
+ (rightPixel
* (unsigned long)(column - *firstRightP)))
/ (unsigned long)v);
} else {
int const v = w * 4;
leftPixel = (sample)(
((leftPixel
* (unsigned long)(Left->pam.width - column))
+ (rightPixel
* (unsigned long)(column - Left->pam.width + v)))
/ (unsigned long)v);
}
}
outPixel[plane] = leftPixel;
}
}
static void
stitchOneRow(Image * const Left,
Image * const Right,
Output * const Out,
Stitcher * const Stitch,
int const row,
int const 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 occurrences
* 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.
*/
tuple * const Row = (*(Out->Row))(Out,row);
unsigned column, firstRight;
firstRight = 0; /* initial value */
for (column = 0; column < Out->image->pam.width; ++column) {
int right_row, right_column;
right_row = -1;
right_column = (*(Stitch->XRight))(Stitch, column,
(y < 0) ? (row + y) : row) + 0.5;
if ((0 <= right_column)
&& (right_column < Right->pam.width)) {
right_row = (*(Stitch->YRight))(Stitch, column,
(y < 0) ? (row + y) : row) + 0.5;
if ((right_row < 0)
|| (Right->pam.height <= right_row)) {
right_column = -1;
right_row = -1;
}
} else
right_column = -1;
/* Create the pixel at column 'column' of row 'row' of the
output 'Out': Row[column].
*/
stitchOnePixel(Left, Right, Out->image->pam, row, column, y,
right_row, right_column, &firstRight, Row[column]);
}
}
static void
stitchit(Image * const Left,
Image * const Right,
const char * const outfilename,
const char * const filter,
Stitcher * const Stitch,
int * const retvalP) {
Output * const Out = findObject(filter, &OutputMethods[0],
sizeof(OutputMethods[0]));
unsigned row;
int xp[8], yp[8], x, y, width, height;
if ((Out == (Output *)NULL) || (Out->Name == (char *)NULL))
*retvalP = -2;
else {
if (verbose)
fprintf (stderr, "Selected %s output filter algorithm\n",
Out->Name);
/* Determine the mask corners for both Left and Right images */
determineMaskCorners(Stitch, Left, Right, xp, yp);
/* Output the combined images */
/* Calculate generic x,y left top corner, and the width and height */
calculateXyWidthHeight(xp, yp, &x, &y, &width, &height);
if (verbose)
printPlan(xp, yp, Left, Right);
if (!(*(Out->Alloc))(Out, outfilename, width, height, &Left->pam))
*retvalP = -9;
else {
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);
}
for (row = 0; row < Out->image->pam.height; row++) {
/* Generate row number 'row' of the output image 'Out' */
stitchOneRow(Left, Right, Out, Stitch, row, y);
(*(Out->FlushRow))(Out,row);
}
(*(Out->FlushImage))(Out);
(*(Out->DeAlloc))(Out);
*retvalP = 0;
}
}
}
static int
pnmstitch(const char * const leftfilename,
const char * const rightfilename,
const char * const outfilename,
int const reqx,
int const reqy,
int const reqWidth,
int const reqHeight,
const char * const stitcher,
const char * const filter)
{
Stitcher * const Stitch = findObject(stitcher, &StitcherMethods[0],
sizeof(StitcherMethods[0]));
Image * Left;
Image * Right;
int retval;
if ((Stitch == (Stitcher *)NULL) || (Stitch->Name == (char *)NULL))
retval = -1;
else {
if (verbose)
fprintf (stderr, "Selected %s stitcher algorithm\n",
Stitch->Name);
/* Left hand image read into memory */
Left = readinit(leftfilename);
if (Left == NULL)
retval = -3;
else {
/* Right hand image read into memory */
Right = readinit(rightfilename);
if (Right == NULL)
retval = -4;
else {
if (Left->pam.depth != Right->pam.depth) {
fprintf(stderr, "Images should have matching depth. "
"The left image has depth %d, "
"while the right has depth %d.",
Left->pam.depth, Right->pam.depth);
retval = -5;
} else if (Left->pam.maxval != Right->pam.maxval) {
fprintf (stderr,
"Images should have matching maxval. "
"The left image has maxval %u, "
"while the right has maxval %u.",
(unsigned)Left->pam.maxval,
(unsigned)Right->pam.maxval);
retval = -6;
} else if ((*(Stitch->Alloc))(Stitch) == FALSE)
retval = -7;
else {
(*(Stitch->Constrain))(Stitch, reqx, reqy,
reqWidth, reqHeight);
if ((*(Stitch->Match))(Stitch, Left, Right) == FALSE)
retval = -8;
else
stitchit(Left, Right, outfilename, filter, Stitch,
&retval);
}
free_image(Right);
}
free_image(Left);
}
}
return retval;
}
/* Output Methods */
/* Helper methods */
static void
OutputDeAlloc(Output * me)
{
if (me->image != (Image *)NULL) {
/* Free up resources */
free_image (me->image);
me->image = (Image *)NULL;
}
if (me->extra != (void *)NULL) {
free (me->extra);
me->extra = (void *)NULL;
}
} /* OutputDeAlloc() - end */
static bool
OutputAlloc(Output * const me,
const char * const file,
unsigned int const width,
unsigned int const height,
struct pam * const 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) {
OutputDeAlloc(me);
return FALSE;
}
return TRUE;
} /* OutputAlloc() - end */
/* StraightThrough output method */
static void
StraightThroughDeAlloc(Output * me)
{
/* Trick the proper freeing of resouces on the Output Image */
me->image->pam.height = 1;
OutputDeAlloc(me);
} /* StraightThroughDeAlloc() - end */
static bool
StraightThroughAlloc(Output * const me,
const char * const file,
unsigned int const width,
unsigned int const height,
struct pam * const prototype)
{
if (OutputAlloc(me, file, width, height, prototype) == FALSE) {
StraightThroughDeAlloc(me);
}
/* Trick the proper allocation of resouces on the Output Image */
me->image->pam.height = 1;
me->image->tuple = pnm_allocpamarray(&me->image->pam);
if (me->image->tuple == (tuple **)NULL) {
StraightThroughDeAlloc(me);
return FALSE;
}
me->image->pam.height = height;
pnm_writepaminit(&me->image->pam);
return TRUE;
} /* StraightThroughAlloc() - end */
static tuple *
StraightThroughRow(Output * me, unsigned row)
{
UNREFERENCED_PARAMETER(row);
return me->image->tuple[0];
} /* StraightThroughRow() - end */
static void
StraightThroughFlushRow(Output * me, unsigned row)
{
UNREFERENCED_PARAMETER(row);
if (me->image != (Image *)NULL) {
pnm_writepamrow(&me->image->pam, me->image->tuple[0]);
}
} /* StraightThroughFlushRow() - end */
static void
StraightThroughFlushImage(Output * me)
{
UNREFERENCED_PARAMETER(me);
} /* StraightThroughFlushImage() - end */
/* Horizontal Crop output method */
#define HorizontalCropDeAlloc StraightThroughDeAlloc
typedef struct {
int state;
int lostInSpace;
} HorizontalCropExtra;
static bool
HorizontalCropAlloc(Output * const me,
const char * const file,
unsigned int const width,
unsigned int const height,
struct pam * const prototype)
{
unsigned long pos;
if (StraightThroughAlloc(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));
/* Test if we can seek, important since we rewrite the header */
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 StraightThroughRow
static void
HorizontalCropFlushRow(Output * me, unsigned row)
{
unsigned column;
unsigned threshold;
# define HorizontalCropThreshold 4
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 = HorizontalCropThreshold;
for (column = 0; column < me->image->pam.width; ++column) {
unsigned plane = 0;
while (me->image->tuple[0][column][plane] == (sample)0) {
if (++plane >= me->image->pam.depth) {
if (--threshold == 0) {
if (((HorizontalCropExtra *)(me->extra))->state == 1) {
((HorizontalCropExtra *)(me->extra))->state = 2;
}
((HorizontalCropExtra *)(me->extra))->lostInSpace++;
return;
}
}
}
if (plane < me->image->pam.depth) {
threshold = HorizontalCropThreshold;
}
}
((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 */
/* Rotate Crop output method */
#define RotateCropDeAlloc OutputDeAlloc
static bool
RotateCropAlloc(Output * const me,
const char * const file,
unsigned int const width,
unsigned int const height,
struct pam * const prototype)
{
if (OutputAlloc(me, file, width, height, prototype) == FALSE) {
RotateCropDeAlloc(me);
}
me->image->tuple = pnm_allocpamarray(&me->image->pam);
if (me->image->tuple == (tuple **)NULL) {
RotateCropDeAlloc(me);
return FALSE;
}
return TRUE;
} /* RotateCropAlloc() - end */
static tuple *
RotateCropRow(Output * me, unsigned row)
{
return me->image->tuple[row];
} /* RotateCropRow() - end */
static void
RotateCropFlushRow(Output * me, unsigned row)
{
UNREFERENCED_PARAMETER(me);
UNREFERENCED_PARAMETER(row);
} /* RotateCropFlushRow() - end */
/*
* Algorithm under construction.
*
*/
static void
RotateCropFlushImage(Output * me)
{
/* Cop Out for now ... */
pnm_writepam(&me->image->pam, me->image->tuple);
} /* RotateCropFlushImage() - end */
/* Output Method Table */
Output OutputMethods[] = {
{ "StraightThrough", StraightThroughAlloc, StraightThroughDeAlloc,
StraightThroughRow, StraightThroughFlushRow,
StraightThroughFlushImage },
{ "HorizontalCrop", HorizontalCropAlloc, HorizontalCropDeAlloc,
HorizontalCropRow, HorizontalCropFlushRow, HorizontalCropFlushImage },
{ "RotateCrop (unimplemented)", RotateCropAlloc, RotateCropDeAlloc,
RotateCropRow, RotateCropFlushRow, RotateCropFlushImage },
{ (char *)NULL }
};
/* Stitcher Methods */
/* These names are for the 8 parameters of a stitch, in any of the 3
methods this program presently implements. Each is a subscript in
the parms[] array for the Stitcher object that represents a linear
stitching method.
There are also other sets of names for the 8 parameters, such as
Rotate_a. I don't know why. Maybe historical.
*/
#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
/* 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)
{
bool retval;
MALLOCARRAY(me->parms, 8);
if (me->parms == NULL)
retval = FALSE;
else {
/* Constraints unset */
me->x = INT_MAX;
me->y = INT_MAX;
me->width = INT_MAX;
me->height = INT_MAX;
/* Unity transform matrix */
me->parms[Sliver_A] = 1.0;
me->parms[Sliver_B] = 0.0;
me->parms[Sliver_C] = 0.0;
me->parms[Sliver_D] = 0.0;
me->parms[Sliver_xp] = 0.0;
me->parms[Sliver_yp] = 1.0;
me->parms[Sliver_xpp] = 0.0;
me->parms[Sliver_ypp] = 0.0;
retval = TRUE;
}
return retval;
}
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 minimum 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,
* but has the possibility of less accuracy.
*
* 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 64
#define SKIP_SLIVER 1
/* Following global variables are for use by SliverMatch() */
static unsigned long starPeriod;
/* The number of events between printing of a "*" progress
indicator.
*/
static unsigned long starCount;
/* The number of events until the next * progress indicator needs to be
printed.
*/
static void
starEvent() {
if (--starCount == 0) {
starCount = starPeriod;
fprintf (stderr, "*");
}
}
static void
starInit(unsigned long const period) {
starPeriod = period;
starCount = period;
}
static void
starResetPeriod(unsigned long const period) {
starPeriod = period;
if (starCount > period)
starCount = period;
}
static void
findBestMatches(Image * const Left,
Image * const Right,
int const x,
int const y,
int const width,
int const height,
int const offY,
unsigned const Xmin,
unsigned const Xmax,
int const Ymin,
int const Ymax,
Best best[NUM_BEST]) {
/*----------------------------------------------------------------------------
Compare the rectangle 'width' columns by 'height' rows with upper
left corner at Column 'x', Row y+offY in image 'Right' to a bunch of
rectangles of the same size in image 'Left' and generate a list of the
rectangles in 'Left' that best match the one in 'Right'.
The specific rectangles in 'Left' we examine are those with upper left
corner (X,Y+offY) where X is in [Xmin, Xmax) and Y is in [Ymin, Ymax).
We return the ordered list of best matches as best[].
Caller must ensure that each of the rectangles in question is fully
contained with its respective image.
-----------------------------------------------------------------------------*/
unsigned X, Y;
/* Exhaustively find the best match */
for (X = Xmin; X < Xmax; ++X) {
int const widthOfOverlap = X - Left->pam.width;
for (Y = Ymin; Y < Ymax; ++Y) {
unsigned long difference = regionDifference(
Left, X, Y + offY,
Right, x, y + offY,
width, height);
update_best(best, difference, widthOfOverlap, Y + offY);
starEvent();
}
}
}
static void
allocate_best_array(Best *** const bestP, unsigned const bestSize) {
Best ** best;
unsigned int i;
MALLOCARRAY(best, bestSize);
if (best == NULL)
pm_error("No memory for Best array");
for (i = 0; i < bestSize; ++i)
best[i] = allocate_best();
*bestP = best;
}
static void determineXYRange(Stitcher * const me,
Image * const Left,
Image * const Right,
unsigned * const XminP,
unsigned * const XmaxP,
int * const YminP,
int * const YmaxP) {
if (me->x == INT_MAX) {
*XmaxP = Left->pam.width - me->width;
/* I can't bring myself to go half way */
*XminP = Left->pam.width - (7 * Right->pam.width / 16);
} else {
*XminP = me->x;
*XmaxP = me->x + 1;
}
if (me->y == INT_MAX) {
/* Middle 1/4 */
*YminP = Left->pam.height * 3 / 8;
*YmaxP = Left->pam.height - (*YminP) - me->height;
} else {
*YminP = me->y;
*YmaxP = me->y + 1;
}
if (verbose)
pm_message("Test %d<x<%d %d<y<%d", *XminP, *XmaxP, *YminP, *YmaxP);
}
/*
* 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 */
unsigned 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;
/* Harry Sticks Geeses */
if (me->width == INT_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);
/* NOTREACHED */
}
if (me->height == INT_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);
/* NOTREACHED */
}
yDiff = (Right->pam.height * skip_sliver) / image_portion;
starInit((unsigned long)-1L);
allocate_best_array(&best, bestSize);
determineXYRange(me, Left, Right, &Xmin, &Xmax, &Ymin, &Ymax);
/* Find the best */
if ((verbose == 1) || (verbose == 2)) {
fprintf (stderr, "%79s|\r|", "");
starInit((unsigned long)
(
(unsigned long)(Xmax - Xmin)
* (unsigned long)(Ymax - Ymin)
) * (unsigned long)bestSize
/ 78L);
}
/* 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.
*/
/*
* First pass through the right hand images to determine
* which are good candidate (top 90 percentile) for content of
* features that we may have a chance of testing with.
*/
{
float minf, maxf;
float * features;
MALLOCARRAY(features, bestSize);
minf = maxf = 0.0;
for (in = 0; in < bestSize; ++in) {
int const offY = yDiff * (in - (bestSize/2));
float SUM[3], SUMSQ[3];
int plane;
for (plane = 0; plane < MIN(Right->pam.depth,3); ++plane) {
SUM[plane] = SUMSQ[plane] = 0.0;
}
for (X = x; X < (x + me->width); ++X) {
for (Y = y + offY; Y < (y + offY + me->height); ++Y) {
for (plane = 0;
plane < MIN(Right->pam.depth,3);
++plane) {
sample point = Right->tuple[Y][X][plane];
SUM[plane] += point;
SUMSQ[plane] += point * point;
}
}
}
/* How many features */
features[in] = 0.0;
for (plane = 0; plane < MIN(Right->pam.depth,3); ++plane) {
features[in] += SUMSQ[plane] -
(SUM[plane]*SUM[plane]/(float)(me->width*me->height));
}
if ((minf == 0.0) || (features[in] < minf)) {
minf = features[in];
}
if ((maxf == 0.0) || (features[in] > maxf)) {
maxf = features[in];
}
}
/* Select 90% in the contrast range */
minf = (minf + maxf) / 10;
for (in = 0; in < bestSize; ++in) {
if (features[in] < minf) {
free_best(best[in]);
best[in] = (Best *)NULL;
}
}
}
/* Loop through the constraints to find the best match */
sumydotx=sumx=sumydot=sumydotydot=sum=sumydoty=sumy=0.0;
xmin = UINT_MAX;
xmax = 0;
ymin = INT_MAX;
ymax = INT_MIN;
in = num = 0;
for (;;) {
float w;
int offY = yDiff * (in - (bestSize/2));
/* See if this one to be skipped because of too few features */
if (best[in] == (Best *)NULL) {
if (in > (bestSize/2)) {
in = bestSize - in;
} else if (in < (bestSize/2)) {
in = (bestSize-1) - in;
} else {
break;
}
if ((verbose == 1) || (verbose == 2))
for (X = Xmin; X < Xmax; ++X) {
for (Y = Ymin; Y < Ymax; ++Y)
starEvent();
}
continue;
}
findBestMatches(Left, Right, x, y, me->width, me->height,
offY, Xmin, Xmax, Ymin, Ymax,
best[in]);
/* slop (noise in NUM_BEST) */
{
float SUMx, SUMxx, SUMy, SUMyy, SUMw;
unsigned i;
SUMx = SUMxx = SUMy = SUMyy = SUMw = 0.0;
for (i = 0; i < NUM_BEST; ++i) {
/* best[in][i] describes the ith closest region in the right
image to the region in the left image whose top corner is
at (x, y+offY).
*/
float const w2 = (best[in][i].total > 0)
? (1.0 / (float)best[in][i].total)
: 1.0;
SUMx += w2 * best[in][i].x;
SUMy += w2 * best[in][i].y;
SUMxx += w2 * (best[in][i].x * best[in][i].x);
SUMyy += w2 * (best[in][i].y * best[in][i].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 > 1) {
if (me->x == INT_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;
}
/* If it is `wacky' drop it */
if ((newXmax - Xmax) < (Right->pam.width / 3)) {
/* 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 == INT_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)) {
starResetPeriod((unsigned long)(
((unsigned long)(Xmax - Xmin)
* (unsigned long)(Ymax - Ymin))
* (unsigned long)bestSize
/ 78L));
}
}
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 (in = 0; in < bestSize; ++in) {
if (best[in] != (Best *)NULL) {
free_best (best[in]);
best[in] = (Best *)NULL;
}
}
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 */
/* These are not used. Perhaps they are forerunners of the more
expressive Sliver_A, etc. names.
*/
#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;
if (verbose)
pm_message("LinearMatch translation parameters are (%d,%d)",
me->x, me->y);
return TRUE;
} /* LinearMatch() - end */
/*
* Transformation parameters
* left x' = x
* left y' = y
* right x' = x + me->x
* right y' = y + me->y
*/
static float
LinearXLeft(Stitcher * me, int x, int y)
{
UNREFERENCED_PARAMETER(y);
return x;
} /* LinearXLeft() - end */
static float
LinearYLeft(Stitcher * me, int x, int y)
{
UNREFERENCED_PARAMETER(x);
return y;
} /* LinearYLeft() - end */
static float
LinearXRight(Stitcher * me, int x, int y)
{
UNREFERENCED_PARAMETER(y);
return (x + me->x);
} /* LinearXRight() - end */
static float
LinearYRight(Stitcher * me, int x, int y)
{
UNREFERENCED_PARAMETER(x);
return (y + me->y);
} /* LinearYRight() - 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 */
static void
BiLinearDeAlloc(Stitcher * me)
{
LinearDeAlloc(me);
}
static bool
BiLinearAlloc(Stitcher * me)
{
return LinearAlloc(me);
}
static void
BiLinearConstrain(Stitcher * me, int x, int y, int width, int height)
{
LinearConstrain(me, x, y, width, height);
if (x != INT_MAX) {
me->parms[3] -= x;
}
if (y != INT_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 minimum 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 xp, yp, xpp, ypp;
if (SliverMatch(me, Left, Right, IMAGE_PORTION, SKIP_SLIVER) == FALSE) {
return FALSE;
}
/* If too wacky, flatten out */
xp = me->parms[Sliver_xp];
yp = me->parms[Sliver_yp];
xpp = me->parms[Sliver_xpp];
ypp = me->parms[Sliver_ypp];
if ((me->parms[Sliver_A] < -0.3)
|| (0.3 < me->parms[Sliver_A])) {
xp = xpp = (xp + xpp) / 2;
}
if ((me->parms[Sliver_C] < 0.6)
|| (1.5 < me->parms[Sliver_D])) {
yp = (yp + ypp - (float)Right->pam.height) / 2;
ypp = yp + Right->pam.height;
}
/*
* 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 float
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];
} /* BiLinearXRight() - end */
static float
BiLinearYRight(Stitcher * me, int x, int y)
{
/* A little trick I learned from a biker */
float 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];
} /* BiLinearYRight() - end */
static void
BiLinearOutput(Stitcher * me, FILE * fp)
{
fprintf (fp,
"x'=%.6gx%+.6gy%+.6gxy%+.6g\ny'=%.6gx'%+.6gy%+.6gx'y%+.6g\n",
me->parms[BiLinear_a], me->parms[BiLinear_b], me->parms[BiLinear_c],
me->parms[BiLinear_d], me->parms[BiLinear_e], me->parms[BiLinear_f],
me->parms[BiLinear_g], me->parms[BiLinear_h]);
} /* BiLinearOutput() - end */
/* Rotate Stitcher Methods */
#define RotateDeAlloc BiLinearDeAlloc
#define RotateAlloc BiLinearAlloc
#define RotateConstrain BiLinearConstrain
/*
* First pass is to utilize the SliverMatch.
*
* Transformation parameters
* x=ax.+by.+d
* y=ex.+fy.+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
*
*/
#define Rotate_a 0
#define Rotate_b 1
#define Rotate_c 2
#define Rotate_d 3
#define Rotate_e 4
#define Rotate_f 5
#define Rotate_g 6
#define Rotate_h 7
static bool
RotateMatch(Stitcher * me, Image * Left, Image * Right)
{
float xp, yp, xpp, ypp;
if (SliverMatch(me, Left, Right, IMAGE_PORTION, SKIP_SLIVER) == FALSE) {
return FALSE;
}
xp = me->parms[Sliver_xp];
yp = me->parms[Sliver_yp];
xpp = me->parms[Sliver_xpp];
ypp = me->parms[Sliver_ypp];
me->parms[Rotate_c] = (xp - xpp);
me->parms[Rotate_c]*= me->parms[Rotate_c];
me->parms[Rotate_g] = (yp - ypp);
me->parms[Rotate_g]*= me->parms[Rotate_g];
me->parms[Rotate_a] = me->parms[Rotate_f] = sqrt(me->parms[Rotate_g]
/ (me->parms[Rotate_c]
+ me->parms[Rotate_g]));
me->parms[Rotate_b] = me->parms[Rotate_e] = sqrt(me->parms[Rotate_c]
/ (me->parms[Rotate_c]
+ me->parms[Rotate_g]));
if (xp < xpp) {
me->parms[Rotate_b] = -me->parms[Rotate_b];
} else {
me->parms[Rotate_e] = -me->parms[Rotate_e];
}
/* negative (for reverse transform below) xp & yp set for unity gain */
xp = ((me->parms[Rotate_b] * (float)Right->pam.height) + xp + xpp) / -2;
yp = ((me->parms[Rotate_f] * (float)Right->pam.height) - yp - ypp) / 2;
me->parms[Rotate_d] = xp * me->parms[Rotate_a] + yp * me->parms[Rotate_b];
me->parms[Rotate_h] = xp * me->parms[Rotate_e] + yp * me->parms[Rotate_f];
return TRUE;
} /* RotateMatch() - end */
/*
* Transformation parameters
* x`=x
* y`=y
*/
#define RotateXLeft BiLinearXLeft
#define RotateYLeft BiLinearYLeft
/*
* Transformation parameters
* x`=ax+by+d
* y`=ex+fy+h
*/
static float
RotateXRight(Stitcher * me, int x, int y)
{
return (me->parms[Rotate_a] * x) + (me->parms[Rotate_b] * y)
+ me->parms[Rotate_d];
} /* RotateXRight() - end */
static float
RotateYRight(Stitcher * me, int x, int y)
{
return (me->parms[Rotate_e] * x) + (me->parms[Rotate_f] * y)
+ me->parms[Rotate_h];
} /* RotateYRight() - end */
static void
RotateOutput(Stitcher * me, FILE * fp)
{
fprintf (fp,
"x'=%.6gx%+.6gy%+.6g\ny'=%.6gx%+.6gy%+.6g\n",
me->parms[Rotate_a], me->parms[Rotate_b], me->parms[Rotate_d],
me->parms[Rotate_e], me->parms[Rotate_f], me->parms[Rotate_h]);
} /* RotateOutput() - end */
/* Stitcher Method Table */
Stitcher StitcherMethods[] = {
{ "RotateSliver", RotateAlloc, RotateDeAlloc, RotateConstrain,
RotateMatch, RotateXLeft, RotateYLeft, RotateXRight,
RotateYRight, RotateOutput },
{ "BiLinearSliver", BiLinearAlloc, BiLinearDeAlloc, BiLinearConstrain,
BiLinearMatch, BiLinearXLeft, BiLinearYLeft, BiLinearXRight,
BiLinearYRight, BiLinearOutput },
{ "LinearSliver", LinearAlloc, LinearDeAlloc, LinearConstrain,
LinearMatch, LinearXLeft, LinearYLeft, LinearXRight,
LinearYRight, LinearOutput },
{ (char *)NULL }
};