{
    if (n < 0) {
	n = strlen (str);
    }

    CRIMP_ASSERT_BOUNDS (n,(buf->sentinel - buf->here));

    if (strncmp(buf->here, (unsigned char*) str, n) != 0) {
    if (strncmp((char*) buf->here, str, n) != 0) {
	return 0;
    }

    buf->here += n;
    return 1;
}

    unsigned char* buf;      /* Start of data */
    unsigned char* here;     /* Current byte, read location */
    unsigned char* sentinel; /* End of buffer, behind last byte */
    int            length;   /* Size of buffer, sentinel - buf */
} crimp_buffer;

#define crimp_buf_at(b) ((b)->here)

/*
 * BUILD ASSERTION: The buffer API assumes that a variable of type 'int' can
 * hold (at least) 4 bytes (See the crimp_read_*int32* functions).  Failure in
 * the line below tells us that this is not true for the chosen combination of
 * OS, compiler, and compiler flags.
 */

CRIMP_BUILD_ASSERT (sizeof(int) >= 4);

/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:

#ifndef CRIMP_COMMON_H
#define CRIMP_COMMON_H
/*
 * CRIMP :: Common Declarations :: PUBLIC
 * (C) 2011.
 * (C) 2011-2012.
 */

/*
 * Support compile time assertions. While mostly intended for type size
 * checks, all C expressions are allowed.
 */
 
#define CRIMP_BUILD_ASSERT(expr)       typedef char CRIMP_BA_UNIQUE_NAME [(expr)?1:-1]
#define CRIMP_BA_UNIQUE_NAME           CRIMP_BA_MAKE_NAME(__LINE__)
#define CRIMP_BA_MAKE_NAME(line)       CRIMP_BA_MAKE_NAME2(line)
#define CRIMP_BA_MAKE_NAME2(line)      __crimp_build_assert_ ## line

/*
 * Checking of 0-based ranges.
 */

#define CRIMP_RANGEOK(i,n) ((0 <= (i)) && (i < (n)))

/*
 * Convenient checking of image types.
 */

    switch(filterPtr->type) {
    case GFT_FIR:
	FIRDestroyFilterSet(&filterPtr->filters.fir);
	break;
    default: 
	break;
    }
    ckfree(filterPtr);
    ckfree((char*)filterPtr);
}

/*
 *-----------------------------------------------------------------------------
 *
 * GaussianFilter01 --
 *

    int width,			/* Width of the images */
    float* inputImage,		/* Input image: (height x width) array of
				 * float's, row-major order */
    float* outputImage		/* Output image: (height x width) array of
				 * float's, row-major order */
) {
    int area = height * width;
    float* tempImage = ckalloc(area * sizeof(float));
    float* tempImage = (float*) ckalloc(area * sizeof(float));

    /* 
     * Filter first the rows and then the columns.
     */

    GaussianFilter01(filterPtr, 0, height, width, inputImage, tempImage);
    GaussianFilter10(filterPtr, 0, height, width, tempImage, outputImage);

    ckfree(tempImage);
    ckfree((char*)tempImage);
}

/*
 *-----------------------------------------------------------------------------
 *
 * GaussianGradientX2D --
 *

    int width,			/* Width of the images */
    float* inputImage,		/* Input image: (height x width) array of
				 * float's, row-major order */
    float* outputImage		/* Output image: (height x width) array of
				 * float's, row-major order */
) {
    int area = height * width;
    float* tempImage = ckalloc(area * sizeof(float));
    float* tempImage = (float*) ckalloc(area * sizeof(float));
    int i;

    /* 
     * Derivative-filter the rows.
     */

    GaussianFilter01(filterPtr, 1, height, width, inputImage, tempImage);

    /*
     * Gaussian-filter the columns
     */

    GaussianFilter10(filterPtr, 0, height, width, tempImage, outputImage);

    ckfree(tempImage);
    ckfree((char*)tempImage);
}

/*
 *-----------------------------------------------------------------------------
 *
 * GaussianGradientY2D --
 *

    int width,			/* Width of the images */
    float* inputImage,		/* Input image: (height x width) array of
				 * float's, row-major order */
    float* outputImage		/* Output image: (height x width) array of
				 * float's, row-major order */
) {
    int area = height * width;
    float* tempImage = ckalloc(area * sizeof(float));
    float* tempImage = (float*) ckalloc(area * sizeof(float));
    int i;

    /* 
     * Gaussian-filter the rows.
     */

    GaussianFilter01(filterPtr, 0, height, width, inputImage, tempImage);

    /*
     * Derivative-filter the columns
     */

    GaussianFilter10(filterPtr, 1, height, width, tempImage, outputImage);

    ckfree(tempImage);
    ckfree((char*)tempImage);
}

/*
 *-----------------------------------------------------------------------------
 *
 * GaussianGradientMagnitude2D --
 *

    float* inputImage,		/* Input image: (height x width) array of
				 * float's, row-major order */
    float* outputImage		/* Output image: (height x width) array of
				 * float's, row-major order */
) {
    int i;
    int area = height * width;
    float* tempImageX = ckalloc(area * sizeof(float));
    float* tempImageX = (float*) ckalloc(area * sizeof(float));

    GaussianGradientX2D(filterPtr, height, width, inputImage, tempImageX);
    GaussianGradientY2D(filterPtr, height, width, inputImage, outputImage);
    for (i = 0; i < area; ++i) {
	outputImage[i] = hypotf(tempImageX[i], outputImage[i]);
    }
    ckfree(tempImageX);
    ckfree((char*)tempImageX);
}

/*
 *-----------------------------------------------------------------------------
 *
 * GaussianLaplacian2D --
 *

    int width,			/* Width of the images */
    float* inputImage,		/* Input image: (height x width) array of
				 * float's, row-major order */
    float* outputImage		/* Output image: (height x width) array of
				 * float's, row-major order */
) {
    int area = height * width;
    float* tempImage1 = ckalloc(area * sizeof(float));
    float* tempImage2 = ckalloc(area * sizeof(float));
    float* tempImage1 = (float*) ckalloc(area * sizeof(float));
    float* tempImage2 = (float*) ckalloc(area * sizeof(float));
    int i;

    /* Gaussian filter by rows */

    GaussianFilter01(filterPtr, 0, height, width, inputImage, tempImage1);

    /* Second-derivative-filter by columns */

    /* Sum the two results */

    for (i = 0; i < area; ++i) {
	outputImage[i] += tempImage2[i];
    }

    ckfree(tempImage1);
    ckfree(tempImage2);
    ckfree((char*)tempImage1);
    ckfree((char*)tempImage2);
}

/*
 *-----------------------------------------------------------------------------
 *
 * FIRInitFilterSet --
 *

    for (i = 0; i < 3; ++i) {
	filterPtr->coefs[i].n0 /= s[i];
	filterPtr->coefs[i].n1 /= s[i];
	filterPtr->coefs[i].n2 /= s[i];
	filterPtr->coefs[i].n3 /= s[i];
    }

    ckfree(y2);
    ckfree(y1);
    ckfree(y0);
    ckfree(x);
    ckfree((char*)y2);
    ckfree((char*)y1);
    ckfree((char*)y0);
    ckfree((char*)x);
}

/*
 *-----------------------------------------------------------------------------
 *
 * DericheApply --
 *

}

crimp_image*
crimp_geo_warp_init (crimp_image* input, crimp_image* forward, int* origx, int* origy)
{
    /*
     * Run the four corners of the input through the forward transformation to
     * get their locations, and use the results to determine dimensions of the
     * output image and the location of its origin point.
     * get their locations, and use the results to determine geometry of the
     * output image, i.e. dimensions and location of its origin point.
     *
     * NOTE: The input image may already come with origin point data. We have
     * to and are taking this into account when computing the input corners.
     * NOTE: We have to take the origin of the input image into account when
     * computing the input corners.
     */

    crimp_image* result;
    double xlu, xru, xld, xrd, left, right;
    double ylu, yru, yld, yrd, up, down;
    int ileft, iright, iup, idown, w, h, iorigx, iorigy, oc = 0;
    Tcl_Obj* meta;
    Tcl_Obj* key1 = Tcl_NewStringObj ("crimp", -1);
    Tcl_Obj* key2 = Tcl_NewStringObj ("origin", -1);
    Tcl_Obj* cmeta = NULL;
    Tcl_Obj* corig = NULL;
    Tcl_Obj* orig [2] = {NULL, NULL};

    if (!input->meta ||	(Tcl_DictObjGet(NULL, input->meta, key1, &cmeta) != TCL_OK) ||
	!cmeta       ||	(Tcl_DictObjGet(NULL, cmeta, key2, &corig)       != TCL_OK) ||
	!corig       || (Tcl_ListObjGetElements(NULL, corig, &oc, &orig) != TCL_OK) ||
	!orig[0]     || (Tcl_GetIntFromObj(NULL,orig[0], &iorigx)        != TCL_OK) ||
    iorigx = crimp_x (input);
	!orig[1]     || (Tcl_GetIntFromObj(NULL,orig[1], &iorigy)        != TCL_OK)) {
	iorigx = iorigy = 0;
    iorigy = crimp_y (input);
    }

    xlu = - iorigx;
    ylu = - iorigy;
    crimp_geo_warp_point (forward, &xlu, &ylu);

    xru = - iorigx + input->w - 1;
    xru = - iorigx + crimp_w(input) - 1;
    yru = - iorigy;
    crimp_geo_warp_point (forward, &xru, &yru);

    xld = - iorigx;
    yld = - iorigy + input->h - 1;
    yld = - iorigy + crimp_h(input);
    crimp_geo_warp_point (forward, &xld, &yld);

    xrd = - iorigx + input->w - 1;
    yrd = - iorigy + input->h - 1;
    xrd = - iorigx + crimp_w(input) - 1;
    yrd = - iorigy + crimp_h(input) - 1;
    crimp_geo_warp_point (forward, &xrd, &yrd);

    left  = MIN (MIN (xlu,xld), MIN (xru,xrd));
    right = MAX (MAX (xlu,xld), MAX (xru,xrd));
    up    = MIN (MIN (ylu,yld), MIN (yru,yrd));
    down  = MAX (MAX (ylu,yld), MAX (yru,yrd));

    ileft  = left;  if (ileft  > left)  ileft --;
    iright = right; if (iright < right) iright ++;
    iup    = up;    if (iup    > up)    iup --;
    idown  = down;  if (idown  < down)  idown ++;

    w = iright - ileft + 1;
    h = idown  - iup   + 1;

    *origx = ileft;
    *origy = iup;

    orig [0] = Tcl_NewIntObj (ileft);
    orig [1] = Tcl_NewIntObj (iup);

    corig = Tcl_NewListObj (2, orig);
    cmeta = Tcl_NewDictObj (); Tcl_DictObjPut (NULL, cmeta, key2, corig);
    meta  = Tcl_NewDictObj (); Tcl_DictObjPut (NULL, meta,  key1, cmeta);
    result = crimp_new_at (input->itype, ileft, iup, w, h);
    return result;
}

extern void
crimp_rect_union (const crimp_geometry* a,
		  const crimp_geometry* b,
		  crimp_geometry* result)
{
    /*
     * Compute the bounding box first, as min and max of the individual
     * boundaries. The max values are one too high, which is canceled
     * when computing the dimensions.
     */

    int minx = MIN (a->x, b->x);
    int miny = MIN (a->y, b->y);
    int maxx = MAX (a->x + a->w, b->x + b->w);
    int maxy = MAX (a->y + a->h, b->y + b->h);

    /*
     * And convert back into a plain geometry with location dimensions.
     */

    return crimp_newm (input->itype, w, h, meta);
    result->x = minx;
    result->y = miny;
    result->w = maxx - minx;
    result->h = maxy - miny;
}

/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

 * API :: Core. 
 */

extern void         crimp_geo_warp_point (crimp_image* matrix, double* x, double* y);
extern crimp_image* crimp_geo_warp_init  (crimp_image* input,
					  crimp_image* forward,
					  int* origx, int* origy);

extern void crimp_rect_union (const crimp_geometry* a,
			      const crimp_geometry* b,
			      crimp_geometry* result);


/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:

};

/*
 * Definitions :: Core.
 */

crimp_image*
crimp_new (const crimp_imagetype* itype, int w, int h)
crimp_new_at (const crimp_imagetype* itype, int x, int y, int w, int h)
{
    /*
     * Note: Pixel storage and header describing it are allocated together.
     */

    size_t       size  = sizeof (crimp_image) + CRIMP_RECT_AREA (w, h) * itype->size;
    crimp_image* image = (crimp_image*) ckalloc (size);

    image->itype = itype;

    image->w     = w;
    image->h     = h;
    image->geo.x = x;
    image->geo.y = y;

    image->geo.w = w;
    image->geo.h = h;

    image->meta  = NULL;

    return image;
}

crimp_image*
crimp_newm (const crimp_imagetype* itype, int w, int h, Tcl_Obj* meta)
crimp_newm_at (const crimp_imagetype* itype, int x, int y, int w, int h, Tcl_Obj* meta)
{
    /*
     * Note: Pixel storage and header describing it are allocated together.
     */

    size_t       size  = sizeof (crimp_image) + CRIMP_RECT_AREA (w, h) * itype->size;
    crimp_image* image = (crimp_image*) ckalloc (size);

    image->itype = itype;

    image->w     = w;
    image->h     = h;
    image->geo.x = x;
    image->geo.y = y;

    image->geo.w = w;
    image->geo.h = h;

    image->meta  = meta;

    if (meta) {
	Tcl_IncrRefCount (meta);
    }

    return image;

    /* image type */
    Tcl_DStringAppendElement (&ds, ci->itype->name);

    /* image width */
    {
	char wstring [20];
	sprintf (wstring, "%u", ci->w);
	sprintf (wstring, "%u", ci->geo.w);
	Tcl_DStringAppendElement (&ds, wstring);
    }

    /* image width */
    {
	char hstring [20];
	sprintf (hstring, "%u", ci->h);
	sprintf (hstring, "%u", ci->geo.h);
	Tcl_DStringAppendElement (&ds, hstring);
    }

    /* image client data */
    if (ci->meta) {
	Tcl_DStringAppendElement (&ds, Tcl_GetString (ci->meta));
    } else {

#ifndef CRIMP_IMAGE_H
#define CRIMP_IMAGE_H
/*
 * CRIMP :: Image Declarations, and API :: PUBLIC
 * (C) 2010 - 2011
 */

#include "common.h"
#include "image_type.h"

/*
 * Structures describing images.
 *
 * - A convenient name for a memory block of pixel data
 * - The geometry (bounding box) of an image.
 * - The image itself.
 */

typedef unsigned char* crimp_pixel_array;

typedef struct crimp_geometry {
    int x; /* Location of the image in the infinite 2D plane */
    int y; /* s.a. */
    int w; /* Image dimension, width  */
    int h; /* Image dimension, height */
} crimp_geometry;

typedef struct crimp_image {
    Tcl_Obj*               meta;     /* Tcl level client data */
    const crimp_imagetype* itype;    /* Reference to type descriptor */
    int                    w;        /* Image dimension, width  */
    crimp_geometry         geo;      /* Image geometry, bounding box */
    int                    h;        /* Image dimension, height */
    unsigned char          pixel[4]; /* Integrated pixel storage */
} crimp_image;

/*
 * Pixel Access Macros. General access to a 'color' channel.
 */

#define CRIMP_CHAN(iptr,c,x,y) ((c) + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define CRIMP_CHAN(iptr,c,x,y) ((c) + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define CH(iptr,c,x,y)   (iptr)->pixel [CRIMP_CHAN (iptr,c,x,y)]

/*
 * Pixel Access Macros. RGBA / RGB
 */

/*
 * Manually optimized, factored the pixelsize out of the summands. It
 * is not sure if this is faster (easier to optimize), or if we should
 * precompute the pitch (w*pixelsize), and have the pixel size mult
 * in each x ... As the pixel size is mostly 1, 2, 4, i.e. redundant
 * removed unity, or a power of 2, i.e handled as shift this should be
 * good enough. The only not so sure case is RGB, with pixel size of 3.
 */

#define SZ(iptr) ((iptr)->itype->size)

#define RED(iptr,x,y)   (0 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define GREEN(iptr,x,y) (1 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define BLUE(iptr,x,y)  (2 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define ALPHA(iptr,x,y) (3 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define RED(iptr,x,y)   (0 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define GREEN(iptr,x,y) (1 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define BLUE(iptr,x,y)  (2 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define ALPHA(iptr,x,y) (3 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))

#if 0 /* Unoptimized formulas */
#define RED(iptr,x,y)   (0 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->w)))
#define GREEN(iptr,x,y) (1 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->w)))
#define BLUE(iptr,x,y)  (2 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->w)))
#define ALPHA(iptr,x,y) (3 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->w)))
#define RED(iptr,x,y)   (0 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->geo.w)))
#define GREEN(iptr,x,y) (1 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->geo.w)))
#define BLUE(iptr,x,y)  (2 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->geo.w)))
#define ALPHA(iptr,x,y) (3 + ((x)*SZ (iptr)) + ((y)*SZ (iptr)*((size_t) (iptr)->geo.w)))
#endif

#define R(iptr,x,y) (iptr)->pixel [RED   (iptr,x,y)]
#define G(iptr,x,y) (iptr)->pixel [GREEN (iptr,x,y)]
#define B(iptr,x,y) (iptr)->pixel [BLUE  (iptr,x,y)]
#define A(iptr,x,y) (iptr)->pixel [ALPHA (iptr,x,y)]

/*
 * Pixel Access Macros. GREY8, GREY16, GREY32, FLOATP.
 *
 * NOTE: The casts should use standard types where we we know the size in
 *       bytes exactly, by definition.
 */

#define CRIMP_INDEX(iptr,x,y) \
    (((x)*SZ (iptr)) + \
     ((y)*SZ (iptr)*(((size_t) (iptr)->w))))
     ((y)*SZ (iptr)*(((size_t) (iptr)->geo.w))))

#define GREY8(iptr,x,y)  *((unsigned char*)  &((iptr)->pixel [CRIMP_INDEX (iptr,x,y)]))
#define GREY16(iptr,x,y) *((unsigned short*) &((iptr)->pixel [CRIMP_INDEX (iptr,x,y)]))
#define GREY32(iptr,x,y) *((unsigned int* )  &((iptr)->pixel [CRIMP_INDEX (iptr,x,y)]))
#define FLOATP(iptr,x,y) *((float*)          &((iptr)->pixel [CRIMP_INDEX (iptr,x,y)]))

/*
 * Pixel as 2-complement numbers (-128..127, instead of unsigned 0..255).
 */

#define SGREY8(iptr,x,y) *((signed char*)  &((iptr)->pixel [CRIMP_INDEX (iptr,x,y)]))

/*
 * Pixel Access Macros. HSV.
 */

#define HUE(iptr,x,y) (0 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define SAT(iptr,x,y) (1 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define VAL(iptr,x,y) (2 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define HUE(iptr,x,y) (0 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define SAT(iptr,x,y) (1 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define VAL(iptr,x,y) (2 + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))

#define H(iptr,x,y) (iptr)->pixel [HUE (iptr,x,y)]
#define S(iptr,x,y) (iptr)->pixel [SAT (iptr,x,y)]
#define V(iptr,x,y) (iptr)->pixel [VAL (iptr,x,y)]

/*
 * Pixel Access Macros. FPCOMPLEX.
 */

#define REAL(iptr,x,y)      (0             + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define IMAGINARY(iptr,x,y) (sizeof(float) + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->w)))
#define REAL(iptr,x,y)      (0             + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))
#define IMAGINARY(iptr,x,y) (sizeof(float) + SZ(iptr) * ((x) + (y)*((size_t) (iptr)->geo.w)))

#define RE(iptr,x,y) *((float*) &((iptr)->pixel [REAL      (iptr,x,y)]))
#define IM(iptr,x,y) *((float*) &((iptr)->pixel [IMAGINARY (iptr,x,y)]))

/*
 * Other constants
 */

#define BLACK 0
#define WHITE 255

#define OPAQUE      255
#define TRANSPARENT 0

/*
 * Area calculations macros.
 */

#define CRIMP_RECT_AREA(w,h) (((size_t) (w)) * (h))
#define crimp_image_area(iptr) (CRIMP_RECT_AREA ((iptr)->w, (iptr)->h))
#define crimp_image_area(iptr) (CRIMP_RECT_AREA (crimp_w(iptr), crimp_h(iptr)))

#define crimp_place(image,ix,iy)			\
    ((image)->geo.x = (ix), (image)->geo.y = (iy))

#define crimp_inside(image,px,py) \
    ((crimp_x(image) <= (px)) && ((px) < (crimp_x(image) + crimp_w(image))) && \
     (crimp_y(image) <= (py)) && ((py) < (crimp_y(image) + crimp_h(image))))

#define crimp_x(image) ((image)->geo.x)
#define crimp_y(image) ((image)->geo.y)
#define crimp_w(image) ((image)->geo.w)
#define crimp_h(image) ((image)->geo.h)

/*
 * Convenience macros for the creation of images with predefined image types.
 */

#define crimp_new_atg(type,g)     (crimp_new_at  ((type), (g).x, (g).y, (g).w, (g).h))
#define crimp_new(type,w,h)       (crimp_new_at  ((type), 0, 0, (w), (h)))
#define crimp_newm(type,w,h,meta) (crimp_newm_at ((type), 0, 0, (w), (h), (meta)))

#define crimp_new_hsv(w,h)       (crimp_new (crimp_imagetype_find ("crimp::image::hsv"),     (w), (h)))
#define crimp_new_rgba(w,h)      (crimp_new (crimp_imagetype_find ("crimp::image::rgba"),    (w), (h)))
#define crimp_new_rgb(w,h)       (crimp_new (crimp_imagetype_find ("crimp::image::rgb"),     (w), (h)))
#define crimp_new_grey8(w,h)     (crimp_new (crimp_imagetype_find ("crimp::image::grey8"),   (w), (h)))
#define crimp_new_grey16(w,h)    (crimp_new (crimp_imagetype_find ("crimp::image::grey16"),  (w), (h)))
#define crimp_new_grey32(w,h)    (crimp_new (crimp_imagetype_find ("crimp::image::grey32"),  (w), (h)))
#define crimp_new_float(w,h)     (crimp_new (crimp_imagetype_find ("crimp::image::float"),   (w), (h)))
#define crimp_new_fpcomplex(w,h) (crimp_new (crimp_imagetype_find ("crimp::image::fpcomplex"), (w), (h)))

#define crimp_new_like(image)           (crimp_newm ((image)->itype, (image)->w, (image)->h, (image)->meta))
#define crimp_new_like_transpose(image) (crimp_newm ((image)->itype, (image)->h, (image)->w, (image)->meta))
#define crimp_new_like(image)           (crimp_newm_at ((image)->itype, crimp_x(image), crimp_y(image), crimp_w(image), crimp_h(image), (image)->meta))
#define crimp_new_like_transpose(image) (crimp_newm_at ((image)->itype, crimp_x(image), crimp_y(image), crimp_h(image), crimp_w(image), (image)->meta))

#define crimp_new_hsv_at(x,y,w,h)       (crimp_new_at (crimp_imagetype_find ("crimp::image::hsv"),       (x), (y), (w), (h)))
#define crimp_new_rgba_at(x,y,w,h)      (crimp_new_at (crimp_imagetype_find ("crimp::image::rgba"),      (x), (y), (w), (h)))
#define crimp_new_rgb_at(x,y,w,h)       (crimp_new_at (crimp_imagetype_find ("crimp::image::rgb"),       (x), (y), (w), (h)))
#define crimp_new_grey8_at(x,y,w,h)     (crimp_new_at (crimp_imagetype_find ("crimp::image::grey8"),     (x), (y), (w), (h)))
#define crimp_new_grey16_at(x,y,w,h)    (crimp_new_at (crimp_imagetype_find ("crimp::image::grey16"),    (x), (y), (w), (h)))
#define crimp_new_grey32_at(x,y,w,h)    (crimp_new_at (crimp_imagetype_find ("crimp::image::grey32"),    (x), (y), (w), (h)))
#define crimp_new_float_at(x,y,w,h)     (crimp_new_at (crimp_imagetype_find ("crimp::image::float"),     (x), (y), (w), (h)))
#define crimp_new_fpcomplex_at(x,y,w,h) (crimp_new_at (crimp_imagetype_find ("crimp::image::fpcomplex"), (x), (y), (w), (h)))

/*
 * Convenience macros for input image handling.
 */

#define crimp_input(objvar,imagevar,itype) \
    if (crimp_get_image_from_obj (interp, (objvar), &(imagevar)) != TCL_OK) { \
	return TCL_ERROR; \
    } \
    CRIMP_ASSERT_IMGTYPE (imagevar, itype)

#define crimp_input_any(objvar,imagevar) \
    if (crimp_get_image_from_obj (interp, (objvar), &(imagevar)) != TCL_OK) { \
	return TCL_ERROR; \
    }

#define crimp_eq_geo(imagea,imageb) \
    (crimp_eq_dim(imagea,imageb) && crimp_eq_loc(imagea,imageb))

#define crimp_eq_dim(imagea,imageb) \
    (crimp_eq_width(imagea,imageb) && crimp_eq_height(imagea,imageb))
    (((imagea)->w == (imageb)->w) && ((imagea)->h == (imageb)->h))

#define crimp_eq_loc(imagea,imageb) \
    (crimp_eq_x(imagea,imageb) && crimp_eq_y(imagea,imageb))

#define crimp_eq_x(imagea,imageb) \
    (crimp_x(imagea) == crimp_x(imageb))

#define crimp_eq_y(imagea,imageb) \
    (crimp_y(imagea) == crimp_y(imageb))

#define crimp_eq_height(imagea,imageb) \
    ((imagea)->h == (imageb)->h)
    (crimp_h(imagea) == crimp_h(imageb))

#define crimp_eq_width(imagea,imageb) \
    ((imagea)->w == (imageb)->w)
    (crimp_w(imagea) == crimp_w(imageb))

#define crimp_require_dim(image,rw,rh)					\
    (((image)->w == (rw)) && ((image)->h == (rh)))
    ((crimp_w(image) == (rw)) && (crimp_h(image) == (rh)))

#define crimp_require_height(image,rh)					\
    ((image)->h == (rh))
    (crimp_h(image) == (rh))

#define crimp_require_width(image,rw)					\
    ((image)->w == (rw))
    (crimp_w(image) == (rw))


/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */
#endif /* CRIMP_IMAGE_H */

    const void* bgValue,
				/* Pointer to the pixel value that will be
				 * used as background. All background pixels
				 * are coalesced into a single component.
				 * NULL means not to use a background value. */
    crimp_image* imagePtr	/* Input image to segment. */
) {
    int height = imagePtr->h;	/* Height of the image */
    int width = imagePtr->w;	/* Width of the image */
    int height = crimp_h(imagePtr);	/* Height of the image */
    int width = crimp_w(imagePtr);	/* Width of the image */
    int locx = crimp_x(imagePtr);	/* Location of the image */
    int locy = crimp_y(imagePtr);	/* Location of the image */
    int esize = SZ(imagePtr);	/* Size of a pixel value */
    int wm1 = width - 1;
    int wp1 = width + 1;
    size_t area = (size_t)width * (size_t)height;
				/* Area of the image in pixels */
    size_t* parent = (size_t*) ckalloc(area * sizeof(size_t));
				/* Parent link data structure for
				 * UNION-FIND partition */
    crimp_image* result = crimp_new_grey32(width, height);
    crimp_image* result = crimp_new_grey32_at (locx, locy, width, height);
				/* Result image containing subset ranks
				 * during the UNION-FIND calculation and
				 * component numbers during the component
				 * numbering pass. */
    RANK_TYPE* rank = (RANK_TYPE*) result->pixel;
				/* Pixel array for the result image */
    int i, j;			/* Row and column indices */

    return result;
}

crimp_image*
crimp_la_multiply_matrix (crimp_image* a, crimp_image* b)
{
    crimp_image* result;
    int x, y, w;
    int x, y, w, n, m;

    CRIMP_ASSERT_IMGTYPE (a, float);
    CRIMP_ASSERT_IMGTYPE (b, float);
    CRIMP_ASSERT (crimp_require_height(a, b->w),"Unable to multiply matrices, size mismatch");
    CRIMP_ASSERT (crimp_require_height(b, a->w),"Unable to multiply matrices, size mismatch");
    CRIMP_ASSERT (crimp_require_height(a, crimp_w(b)),"Unable to multiply matrices, size mismatch");
    CRIMP_ASSERT (crimp_require_height(b, crimp_w(a)),"Unable to multiply matrices, size mismatch");

    n = crimp_h (a);
    m = crimp_w (a);

    result = crimp_new_float (a->h, a->h);
    result = crimp_new_float (n, n);

    for (y = 0; y < a->h; y++) {
	for (x = 0; x < a->h; x++) {
    for (y = 0; y < n; y++) {
	for (x = 0; x < n; x++) {

	    FLOATP (result, x, y) = 0;
	    for (w = 0; w < a->w; w++) {
	    for (w = 0; w < m; w++) {
		FLOATP (result, x, y) += FLOATP (a, w, y) * FLOATP (b, x, w);
	    }
	}
    }

    return result;
}

};

/*
 * Definitions :: Core.
 */

crimp_volume*
crimp_vnew (const crimp_imagetype* itype, int w, int h, int d)
crimp_vnew_at (const crimp_imagetype* itype, int x, int y, int z, int w, int h, int d)
{
    /*
     * Note: Pixel storage and header describing it are allocated together.
     */

    size_t        size   = sizeof (crimp_volume) + CRIMP_RECT_VOLUME (w, h, d) * itype->size;
    crimp_volume* volume = (crimp_volume*) ckalloc (size);

    volume->itype = itype;

    volume->w     = w;
    volume->h     = h;
    volume->d     = d;
    volume->geo.x = x;
    volume->geo.y = y;
    volume->geo.z = z;

    volume->geo.w = w;
    volume->geo.h = h;
    volume->geo.d = d;

    volume->meta  = NULL;

    return volume;
}

crimp_volume*
crimp_vnewm (const crimp_imagetype* itype, int w, int h, int d, Tcl_Obj* meta)
crimp_vnewm_at (const crimp_imagetype* itype, int x, int y, int z, int w, int h, int d, Tcl_Obj* meta)
{
    /*
     * Note: Pixel storage and header describing it are allocated together.
     */

    size_t        size   = sizeof (crimp_volume) + CRIMP_RECT_VOLUME (w, h, d) * itype->size;
    crimp_volume* volume = (crimp_volume*) ckalloc (size);

    volume->itype = itype;

    volume->w     = w;
    volume->h     = h;
    volume->d     = d;
    volume->geo.x = x;
    volume->geo.y = y;
    volume->geo.z = z;

    volume->geo.w = w;
    volume->geo.h = h;
    volume->geo.d = d;

    volume->meta  = meta;

    if (meta) {
	Tcl_IncrRefCount (meta);
    }

    return volume;

    /* volume type */
    Tcl_DStringAppendElement (&ds, cv->itype->name);

    /* volume width */
    {
	char wstring [20];
	sprintf (wstring, "%u", cv->w);
	sprintf (wstring, "%u", cv->geo.w);
	Tcl_DStringAppendElement (&ds, wstring);
    }

    /* volume height */
    {
	char hstring [20];
	sprintf (hstring, "%u", cv->h);
	sprintf (hstring, "%u", cv->geo.h);
	Tcl_DStringAppendElement (&ds, hstring);
    }

    /* volume depth */
    {
	char dstring [20];
	sprintf (dstring, "%u", cv->d);
	sprintf (dstring, "%u", cv->geo.d);
	Tcl_DStringAppendElement (&ds, dstring);
    }

    /* volume client data */
    if (cv->meta) {
	Tcl_DStringAppendElement (&ds, Tcl_GetString (cv->meta));
    } else {

#ifndef CRIMP_VOLUME_H
#define CRIMP_VOLUME_H
/*
 * CRIMP :: Volume Declarations, and API :: PUBLIC
 * (C) 2010 - 2011
 */

#include "common.h"
#include "image_type.h"

/*
 * Structures describing volumes.
 *
 * - The geometry (bounding box) of a volume.
 * - The volume itself.
 */

typedef struct crimp_geometry3d {
    int x; /* Location of the volume in the infinite 3D volume */
    int y; /* s.a. */
    int z; /* s.a. */
    int w; /* Volume dimension, width  */
    int h; /* Volume dimension, height */
    int d; /* Volume dimension, depth */
} crimp_geometry3d;

typedef struct crimp_volume {
    Tcl_Obj*               meta;     /* Tcl level client data */
    const crimp_imagetype* itype;    /* Reference to type descriptor */
    int                    w;        /* Volume dimension, width  */
    crimp_geometry3d       geo;      /* Volume geometry, bounding box */
    int                    h;        /* Volume dimension, height */
    int                    d;        /* Volume dimension, depth */
    unsigned char          voxel[4]; /* Integrated voxel storage */
} crimp_volume;

/*
 * Voxel Access Macros.
 */

#define CRIMP_VINDEX(iptr,x,y,z) \
    (((x)*SZ (iptr)) + \
     ((y)*SZ (iptr)*((iptr)->w)) + \
     ((z)*SZ (iptr)*((iptr)->w)*((size_t) (iptr)->h)))
     ((y)*SZ (iptr)*((iptr)->geo.w)) + \
     ((z)*SZ (iptr)*((iptr)->geo.w)*((size_t) (iptr)->geo.h)))

#define VFLOATP(iptr,x,y,z) *((float*) &((iptr)->voxel [CRIMP_VINDEX (iptr,x,y,z)]))

/*
 * Convenience macros for the creation of volumes with predefined image types.
 */

#define crimp_vnew_atg(type,g)       (crimp_vnew_at  ((type), (g).x, (g).y, (g).z, (g).w, (g).h, (g).d))
#define crimp_vnew(type,w,h,d)       (crimp_vnew_at  ((type), 0, 0, 0, (w), (h), (d)))
#define crimp_vnewm(type,w,h,d,meta) (crimp_vnewm_at ((type), 0, 0, 0, (w), (h), (d), (meta)))

#define crimp_vnew_hsv(w,h,d)       (crimp_vnew (crimp_imagetype_find ("crimp::image::hsv"),     (w), (h), (d)))
#define crimp_vnew_rgba(w,h,d)      (crimp_vnew (crimp_imagetype_find ("crimp::image::rgba"),    (w), (h), (d)))
#define crimp_vnew_rgb(w,h,d)       (crimp_vnew (crimp_imagetype_find ("crimp::image::rgb"),     (w), (h), (d)))
#define crimp_vnew_grey8(w,h,d)     (crimp_vnew (crimp_imagetype_find ("crimp::image::grey8"),   (w), (h), (d)))
#define crimp_vnew_grey16(w,h,d)    (crimp_vnew (crimp_imagetype_find ("crimp::image::grey16"),  (w), (h), (d)))
#define crimp_vnew_grey32(w,h,d)    (crimp_vnew (crimp_imagetype_find ("crimp::image::grey32"),  (w), (h), (d)))
#define crimp_vnew_float(w,h,d)     (crimp_vnew (crimp_imagetype_find ("crimp::image::float"),   (w), (h), (d)))
#define crimp_vnew_fpcomplex(w,h,d) (crimp_vnew (crimp_imagetype_find ("crimp::image::fpcomplex"), (w), (h), (d)))

#define crimp_vnew_like(volume)           (crimp_vnewm ((volume)->itype, (volume)->w, (volume)->h, (volume)->d, (volume)->meta))
#define crimp_vnew_like_transpose(volume) (crimp_vnewm ((volume)->itype, (volume)->h, (volume)->w, (volume)->d, (volume)->meta))
#define crimp_vnew_like(volume)           (crimp_vnewm ((volume)->itype, (volume)->geo.w, (volume)->geo.h, (volume)->geo.d, (volume)->meta))
#define crimp_vnew_like_transpose(volume) (crimp_vnewm ((volume)->itype, (volume)->geo.h, (volume)->geo.w, (volume)->geo.d, (volume)->meta))

/*
 * Volume calculations macros.
 */

#define CRIMP_RECT_VOLUME(w,h,d) (((size_t) (w)) * (h) * (d))
#define crimp_volume_vol(vptr) (CRIMP_RECT_VOLUME ((vptr)->w, (vptr)->h, (vptr)->d))
#define crimp_volume_vol(vptr) (CRIMP_RECT_VOLUME ((vptr)->geo.w, (vptr)->geo.h, (vptr)->geo.d))

#define crimp_vplace(image,ix,iy,iz)			\
    ((image)->geo.x = (ix), (image)->geo.y = (iy), (image)->geo.z = (iz))

#define crimp_vinside(volume,px,py,pz)					\
    (((volume)->geo.x <= (px)) && ((px) < ((volume)->geo.x + (volume)->geo.w)) && \
     ((volume)->geo.y <= (py)) && ((py) < ((volume)->geo.y + (volume)->geo.h)) && \
     ((volume)->geo.z <= (pz)) && ((pz) < ((volume)->geo.z + (volume)->geo.d)))

#define crimp_vx(volume) ((volume)->geo.x)
#define crimp_vy(volume) ((volume)->geo.y)
#define crimp_vz(volume) ((volume)->geo.z)
#define crimp_vw(volume) ((volume)->geo.w)
#define crimp_vh(volume) ((volume)->geo.h)
#define crimp_vd(volume) ((volume)->geo.d)

/*
 * Convenience macros for input volume handling.
 */

#define crimp_vinput(objvar,volumevar,itype) \
    if (crimp_get_volume_from_obj (interp, (objvar), &(volumevar)) != TCL_OK) { \

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, float);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	FLOATP (result, x, y) = BINOP (FLOATP (imageA, x, y), FLOATP (imageB, x, y));
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, float);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_v = ina ? FLOATP (imageA, pxa, pya) : BLACK;
	double b_v = inb ? FLOATP (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, fpcomplex);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_v = ina ? FLOATP (imageA, pxa, pya) : BLACK;
	double b_re = inb ? RE (imageB, pxb, pyb) : BLACK;
	double b_im = inb ? IM (imageB, pxb, pyb) : BLACK;
	
	double z_vre = BINOP (a_v, b_re);
	double z_vim = BINOP (a_v, b_im);
	
	RE (result, px, py) = BINOP_POST (z_vre);
	IM (result, px, py) = BINOP_POST (z_vim);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, grey16);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	FLOATP (result, x, y) = BINOP (FLOATP (imageA, x, y), GREY16 (imageB, x, y));
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, grey16);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_v = ina ? FLOATP (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY16 (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, grey32);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	FLOATP (result, x, y) = BINOP (FLOATP (imageA, x, y), GREY32 (imageB, x, y));
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, grey32);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_v = ina ? FLOATP (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY32 (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, grey8);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	FLOATP (result, x, y) = BINOP (FLOATP (imageA, x, y), GREY8 (imageB, x, y));
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, float);
crimp_input (imageBObj, imageB, grey8);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_v = ina ? FLOATP (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, fpcomplex);
crimp_input (imageBObj, imageB, float);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_re = ina ? RE (imageA, pxa, pya) : BLACK;
	double a_im = ina ? IM (imageA, pxa, pya) : BLACK;
	double b_v = inb ? FLOATP (imageB, pxb, pyb) : BLACK;
	
	double z_rev = BINOP (a_re, b_v);
	double z_imv = BINOP (a_im, b_v);
	
	RE (result, px, py) = BINOP_POST (z_rev);
	IM (result, px, py) = BINOP_POST (z_imv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA,fpcomplex);
crimp_input (imageBObj, imageB,fpcomplex);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	RE (result, x, y) = BINOP (RE (imageA, x, y), RE (imageB, x, y));
	IM (result, x, y) = BINOP (IM (imageA, x, y), IM (imageB, x, y));
	}
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, fpcomplex);
crimp_input (imageBObj, imageB, fpcomplex);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_re = ina ? RE (imageA, pxa, pya) : BLACK;
	double a_im = ina ? IM (imageA, pxa, pya) : BLACK;
	double b_re = inb ? RE (imageB, pxb, pyb) : BLACK;
	double b_im = inb ? IM (imageB, pxb, pyb) : BLACK;
	
	double z_rere = BINOP (a_re, b_re);
	double z_imim = BINOP (a_im, b_im);
	
	RE (result, px, py) = BINOP_POST (z_rere);
	IM (result, px, py) = BINOP_POST (z_imim);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, fpcomplex);
crimp_input (imageBObj, imageB, fpcomplex);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_GLOBAL
#define BINOP_GLOBAL(a,b,c,d)
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_re = ina ? RE (imageA, pxa, pya) : BLACK;
	double a_im = ina ? IM (imageA, pxa, pya) : BLACK;
	double b_re = inb ? RE (imageB, pxb, pyb) : BLACK;
	double b_im = inb ? IM (imageB, pxb, pyb) : BLACK;

	BINOP_GLOBAL (a_re, a_im, b_re, b_im);
	
	RE (result, px, py) = BINOP_RE (a_re, a_im, b_re, b_im);
	IM (result, px, py) = BINOP_IM (a_re, a_im, b_re, b_im);
    }
}

#undef BINOP_RE
#undef BINOP_IM
#undef BINOP_GLOBAL

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, fpcomplex);
crimp_input (imageBObj, imageB, grey16);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_re = ina ? RE (imageA, pxa, pya) : BLACK;
	double a_im = ina ? IM (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY16 (imageB, pxb, pyb) : BLACK;
	
	double z_rev = BINOP (a_re, b_v);
	double z_imv = BINOP (a_im, b_v);
	
	RE (result, px, py) = BINOP_POST (z_rev);
	IM (result, px, py) = BINOP_POST (z_imv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, fpcomplex);
crimp_input (imageBObj, imageB, grey32);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_re = ina ? RE (imageA, pxa, pya) : BLACK;
	double a_im = ina ? IM (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY32 (imageB, pxb, pyb) : BLACK;
	
	double z_rev = BINOP (a_re, b_v);
	double z_imv = BINOP (a_im, b_v);
	
	RE (result, px, py) = BINOP_POST (z_rev);
	IM (result, px, py) = BINOP_POST (z_imv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, fpcomplex);
crimp_input (imageBObj, imageB, grey8);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	double a_re = ina ? RE (imageA, pxa, pya) : BLACK;
	double a_im = ina ? IM (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (imageB, pxb, pyb) : BLACK;
	
	double z_rev = BINOP (a_re, b_v);
	double z_imv = BINOP (a_im, b_v);
	
	RE (result, px, py) = BINOP_POST (z_rev);
	IM (result, px, py) = BINOP_POST (z_imv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, float);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	FLOATP (result, x, y) = BINOP (GREY16 (imageA, x, y), FLOATP (imageB, x, y));
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, float);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	double b_v = inb ? FLOATP (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, float);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_grey16_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	double b_v = inb ? FLOATP (imageB, pxb, pyb) : BLACK;
	
	int    z_vv = BINOP (a_v, b_v);
	
	GREY16 (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, fpcomplex);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_fpcomplex_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	double b_re = inb ? RE (imageB, pxb, pyb) : BLACK;
	double b_im = inb ? IM (imageB, pxb, pyb) : BLACK;
	
	double z_vre = BINOP (a_v, b_re);
	double z_vim = BINOP (a_v, b_im);
	
	RE (result, px, py) = BINOP_POST (z_vre);
	IM (result, px, py) = BINOP_POST (z_vim);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int x, y;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, grey16);

if (!crimp_eq_dim (imageA, imageB)) {
    Tcl_SetResult(interp, "image dimensions do not match", TCL_STATIC);
    return TCL_ERROR;
}

result = crimp_new_like (imageA);

for (y = 0; y < result->h; y++) {
    for (x = 0; x < result->w; x++) {

	GREY16 (result, x, y) = BINOP (GREY16 (imageA, x, y), GREY16 (imageB, x, y));
    }
}

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, grey16);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY16 (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, grey16);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_grey16_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY16 (imageB, pxb, pyb) : BLACK;
	
	int    z_vv = BINOP (a_v, b_v);
	
	GREY16 (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, grey32);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY32 (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, grey32);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_grey32_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY32 (imageB, pxb, pyb) : BLACK;
	
	int    z_vv = BINOP (a_v, b_v);
	
	GREY32 (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */

crimp_image*     result;
crimp_image*     imageA;
crimp_image*     imageB;
int px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey16);
crimp_input (imageBObj, imageB, grey8);

/*
 * Compute union area of the two images to process.
 * Note how the images do not have to match in size.
 */

crimp_rect_union (&imageA->geo, &imageB->geo, &bb);

result = crimp_new_float_at (bb.x, bb.y, bb.w, bb.h);
oxa = crimp_x (imageA);
oya = crimp_y (imageA);
oxb = crimp_x (imageB);
oyb = crimp_y (imageB);

/*
 * px, py are physical coordinates in the result, starting from 0.
 * The associated logical coordinates in the 2D plane are
 *  lx = px + x(result)
 *  lx = py + y(result)
 * And when we are inside an input its physical coordinates, from the logical are
 *  px' = lx - x(input)
 *  py' = ly - y(input)
 * Important to note, we can compute all these directly as loop variables, as
 * they are all linearly related to each other.
 */

#ifndef BINOP_POST
#define BINOP_POST(z) z
#endif

for (py = 0, ly = bb.y, pya = bb.y - oya, pyb = bb.y - oyb;
     py < bb.h;
     py++, ly++, pya++, pyb++) {

    for (px = 0, lx = bb.x, pxa = bb.x - oxa, pxb = bb.x - oxb;
	 px < bb.w;
	 px++, lx++, pxa++, pxb++) {

	int ina = crimp_inside (imageA, lx, ly);
	int inb = crimp_inside (imageB, lx, ly);

	/*
	 * The result depends on where we are relative to both input.
	 * Inside of each input we take the respective value of the
	 * pixel. Outside of an input we take BLACK as the value
	 * instead.
	 */

	int    a_v = ina ? GREY16 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (imageB, pxb, pyb) : BLACK;
	
	double z_vv = BINOP (a_v, b_v);
	
	FLOATP (result, px, py) = BINOP_POST (z_vv);
    }
}

#undef BINOP
#undef BINOP_POST

Tcl_SetObjResult(interp, crimp_new_image_obj (result));
return TCL_OK;

/* vim: set sts=4 sw=4 tw=80 et ft=c: */
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */