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

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

    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-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((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 = (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((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 = (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((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 = (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((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 = (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((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 = (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((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((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 geometry of the
     * output image, i.e. dimensions and location of its origin point.
     *
     * 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;










    iorigx = crimp_x (input);

    iorigy = crimp_y (input);


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

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

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

    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;

    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.
     */

    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_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->geo.x = x;
    image->geo.y = y;

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

    image->meta  = NULL;

    return image;
}

crimp_image*
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->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->geo.w);
	Tcl_DStringAppendElement (&ds, wstring);
    }

    /* image width */
    {
	char hstring [20];
	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 */
    crimp_geometry         geo;      /* Image geometry, bounding box */

    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)->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)->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)->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)->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)->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)->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 (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_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))

#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) \
    (crimp_h(imagea) == crimp_h(imageb))

#define crimp_eq_width(imagea,imageb) \
    (crimp_w(imagea) == crimp_w(imageb))

#define crimp_require_dim(image,rw,rh)					\
    ((crimp_w(image) == (rw)) && (crimp_h(image) == (rh)))

#define crimp_require_height(image,rh)					\
    (crimp_h(image) == (rh))

#define crimp_require_width(image,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 = 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_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, n, m;

    CRIMP_ASSERT_IMGTYPE (a, float);
    CRIMP_ASSERT_IMGTYPE (b, float);
    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 (n, n);

    for (y = 0; y < n; y++) {
	for (x = 0; x < n; x++) {

	    FLOATP (result, x, y) = 0;
	    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_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->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_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->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->geo.w);
	Tcl_DStringAppendElement (&ds, wstring);
    }

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

    /* volume depth */
    {
	char dstring [20];
	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 */
    crimp_geometry3d       geo;      /* Volume geometry, bounding box */


    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)->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)->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)->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 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 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 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 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 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 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 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:
 */

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_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 ? GREY8 (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, grey32);
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 ? GREY32 (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, grey32);
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_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 ? GREY32 (imageA, pxa, pya) : BLACK;
	double b_v = inb ? FLOATP (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, grey32);
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 ? GREY32 (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 px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey32);
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 ? GREY32 (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, grey32);
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_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 ? GREY32 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY16 (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, grey32);
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 ? GREY32 (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, grey32);
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 ? GREY32 (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, grey32);
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 ? GREY32 (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, grey32);
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_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 ? GREY32 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (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, grey8);
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 ? GREY8 (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, grey8);
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_grey8_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 ? GREY8 (imageA, pxa, pya) : BLACK;
	double b_v = inb ? FLOATP (imageB, pxb, pyb) : BLACK;
	
	int    z_vv = BINOP (a_v, b_v);
	
	GREY8 (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, grey8);
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 ? GREY8 (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 px, py, lx, ly, oxa, oya, oxb, oyb;
int pxa, pya, pxb, pyb;
crimp_geometry bb;

crimp_input (imageAObj, imageA, grey8);
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 ? GREY8 (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, grey8);
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 ? GREY8 (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, grey8);
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 ? GREY8 (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, grey8);
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 ? GREY8 (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, grey8);
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 ? GREY8 (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, grey8);
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_grey8_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 ? GREY8 (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (imageB, pxb, pyb) : BLACK;
	
	int    z_vv = BINOP (a_v, b_v);
	
	GREY8 (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, grey8);
crimp_input (imageBObj, imageB, rgb);

/*
 * 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_rgb_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 ? GREY8 (imageA, pxa, pya) : BLACK;
	int    b_r = inb ? R (imageB, pxb, pyb) : BLACK;
	int    b_g = inb ? G (imageB, pxb, pyb) : BLACK;
	int    b_b = inb ? B (imageB, pxb, pyb) : BLACK;
	
	int    z_vr = BINOP (a_v, b_r);
	int    z_vg = BINOP (a_v, b_g);
	int    z_vb = BINOP (a_v, b_b);
	
	R (result, px, py) = BINOP_POST (z_vr);
	G (result, px, py) = BINOP_POST (z_vg);
	B (result, px, py) = BINOP_POST (z_vb);
    }
}

#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, grey8);
crimp_input (imageBObj, imageB, rgba);

/*
 * 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_rgba_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 ? GREY8 (imageA, pxa, pya) : BLACK;
	int    a_a = ina ? OPAQUE : BLACK;
	int    b_r = inb ? R (imageB, pxb, pyb) : BLACK;
	int    b_g = inb ? G (imageB, pxb, pyb) : BLACK;
	int    b_b = inb ? B (imageB, pxb, pyb) : BLACK;
	int    b_a = inb ? A (imageB, pxb, pyb) : BLACK;
	
	R (result, px, py) = BINOP (a_v, b_r);
	G (result, px, py) = BINOP (a_v, b_g);
	B (result, px, py) = BINOP (a_v, b_b);
	A (result, px, py) = MAX   (a_a, b_a);
    }
}

#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, hsv);
crimp_input (imageBObj, imageB, hsv);

/*
 * 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_hsv_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_h = ina ? H (imageA, pxa, pya) : BLACK;
	int    a_s = ina ? S (imageA, pxa, pya) : BLACK;
	int    a_v = ina ? V (imageA, pxa, pya) : BLACK;
	int    b_h = inb ? H (imageB, pxb, pyb) : BLACK;
	int    b_s = inb ? S (imageB, pxb, pyb) : BLACK;
	int    b_v = inb ? V (imageB, pxb, pyb) : BLACK;
	
	int    z_hh = BINOP (a_h, b_h);
	int    z_ss = BINOP (a_s, b_s);
	int    z_vv = BINOP (a_v, b_v);
	
	H (result, px, py) = BINOP_POST (z_hh);
	S (result, px, py) = BINOP_POST (z_ss);
	V (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, rgb);
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_rgb_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_r = ina ? R (imageA, pxa, pya) : BLACK;
	int    a_g = ina ? G (imageA, pxa, pya) : BLACK;
	int    a_b = ina ? B (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (imageB, pxb, pyb) : BLACK;
	
	int    z_rv = BINOP (a_r, b_v);
	int    z_gv = BINOP (a_g, b_v);
	int    z_bv = BINOP (a_b, b_v);
	
	R (result, px, py) = BINOP_POST (z_rv);
	G (result, px, py) = BINOP_POST (z_gv);
	B (result, px, py) = BINOP_POST (z_bv);
    }
}

#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, rgb);
crimp_input (imageBObj, imageB, rgb);

/*
 * 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_rgb_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_r = ina ? R (imageA, pxa, pya) : BLACK;
	int    a_g = ina ? G (imageA, pxa, pya) : BLACK;
	int    a_b = ina ? B (imageA, pxa, pya) : BLACK;
	int    b_r = inb ? R (imageB, pxb, pyb) : BLACK;
	int    b_g = inb ? G (imageB, pxb, pyb) : BLACK;
	int    b_b = inb ? B (imageB, pxb, pyb) : BLACK;
	
	int    z_rr = BINOP (a_r, b_r);
	int    z_gg = BINOP (a_g, b_g);
	int    z_bb = BINOP (a_b, b_b);
	
	R (result, px, py) = BINOP_POST (z_rr);
	G (result, px, py) = BINOP_POST (z_gg);
	B (result, px, py) = BINOP_POST (z_bb);
    }
}

#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, rgb);
crimp_input (imageBObj, imageB, rgba);

/*
 * 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_rgba_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_r = ina ? R (imageA, pxa, pya) : BLACK;
	int    a_g = ina ? G (imageA, pxa, pya) : BLACK;
	int    a_b = ina ? B (imageA, pxa, pya) : BLACK;
	int    a_a = ina ? OPAQUE : BLACK;
	int    b_r = inb ? R (imageB, pxb, pyb) : BLACK;
	int    b_g = inb ? G (imageB, pxb, pyb) : BLACK;
	int    b_b = inb ? B (imageB, pxb, pyb) : BLACK;
	int    b_a = inb ? A (imageB, pxb, pyb) : BLACK;
	
	R (result, px, py) = BINOP (a_r, b_r);
	G (result, px, py) = BINOP (a_g, b_g);
	B (result, px, py) = BINOP (a_b, b_b);
	A (result, px, py) = MAX   (a_a, b_a);
    }
}

#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, rgba);
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_rgba_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_r = ina ? R (imageA, pxa, pya) : BLACK;
	int    a_g = ina ? G (imageA, pxa, pya) : BLACK;
	int    a_b = ina ? B (imageA, pxa, pya) : BLACK;
	int    a_a = ina ? A (imageA, pxa, pya) : BLACK;
	int    b_v = inb ? GREY8 (imageB, pxb, pyb) : BLACK;
	int    b_a = inb ? OPAQUE : BLACK;
	
	R (result, px, py) = BINOP (a_r, b_v);
	G (result, px, py) = BINOP (a_g, b_v);
	B (result, px, py) = BINOP (a_b, b_v);
	A (result, px, py) = MAX   (a_a, b_a);
    }
}

#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, rgba);
crimp_input (imageBObj, imageB, rgb);

/*
 * 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_rgba_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_r = ina ? R (imageA, pxa, pya) : BLACK;
	int    a_g = ina ? G (imageA, pxa, pya) : BLACK;
	int    a_b = ina ? B (imageA, pxa, pya) : BLACK;
	int    a_a = ina ? A (imageA, pxa, pya) : BLACK;
	int    b_r = inb ? R (imageB, pxb, pyb) : BLACK;
	int    b_g = inb ? G (imageB, pxb, pyb) : BLACK;
	int    b_b = inb ? B (imageB, pxb, pyb) : BLACK;
	int    b_a = inb ? OPAQUE : BLACK;
	
	R (result, px, py) = BINOP (a_r, b_r);
	G (result, px, py) = BINOP (a_g, b_g);
	B (result, px, py) = BINOP (a_b, b_b);
	A (result, px, py) = MAX   (a_a, b_a);
    }
}

#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, rgba);
crimp_input (imageBObj, imageB, rgba);

/*
 * 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_rgba_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_r = ina ? R (imageA, pxa, pya) : BLACK;
	int    a_g = ina ? G (imageA, pxa, pya) : BLACK;
	int    a_b = ina ? B (imageA, pxa, pya) : BLACK;
	int    a_a = ina ? A (imageA, pxa, pya) : BLACK;
	int    b_r = inb ? R (imageB, pxb, pyb) : BLACK;
	int    b_g = inb ? G (imageB, pxb, pyb) : BLACK;
	int    b_b = inb ? B (imageB, pxb, pyb) : BLACK;
	int    b_a = inb ? A (imageB, pxb, pyb) : BLACK;
	
	R (result, px, py) = BINOP (a_r, b_r);
	G (result, px, py) = BINOP (a_g, b_g);
	B (result, px, py) = BINOP (a_b, b_b);
	A (result, px, py) = MAX   (a_a, b_a);
    }
}

#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, @TYPE_A@);
crimp_input (imageBObj, imageB, @TYPE_B@);

/*
 * 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_@TYPE_Z@_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.
	 */

@TRANSFORM@
    }
}

#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:
 */

    int          xo, yo, xi, yi;

    if ((ww < 0) || (hn < 0) || (we < 0) || (hs < 0)) {
	Tcl_SetResult(interp, "bad image border size, expected non-negative values", TCL_STATIC);
	return TCL_ERROR;
    }

    result = crimp_new_at (image->itype, 
			   crimp_x (image) - ww,
			   crimp_y (image) - hn,
			   crimp_w (image) + ww + we,
			   crimp_h (image) + hn + hs);

    /*
     * Nine quadrants to fill:
     *
     * NW N NE
     *  W C E
     * SW S SE

    /*
     * North.
     */

    if (hn) {
	for (yo = 0; yo < hn; yo++) {
	    for (xo = 0; xo < crimp_w (image); xo++) {
		FILL_N (xo + ww, yo);
	    }
	}
    }

    /*
     * North East.
     */

    if (hn && we) {
	for (yo = 0; yo < hn; yo++) {
	    for (xo = 0; xo < we; xo++) {
		FILL_NE (xo + crimp_w (image) + ww, yo);
	    }
	}
    }

    /*
     * West.
     */

    if (ww) {
	for (xo = 0; xo < ww; xo++) {
	    for (yo = 0; yo < crimp_h (image); yo++) {
		FILL_W (xo, yo + hn);
	    }
	}
    }

    /*
     * East.
     */

    if (we) {
	for (xo = 0; xo < we; xo++) {
	    for (yo = 0; yo < crimp_h (image); yo++) {
		FILL_E (xo + crimp_w (image) + ww, yo + hn);
	    }
	}
    }

    /*
     * South West.
     */

    if (hs && ww) {
	for (yo = 0; yo < hs; yo++) {
	    for (xo = 0; xo < ww; xo++) {
		FILL_SW (xo, yo + crimp_h (image) + hn);
	    }
	}
    }

    /*
     * South.
     */

    if (hs) {
	for (yo = 0; yo < hs; yo++) {
	    for (xo = 0; xo < crimp_w (image); xo++) {
		FILL_S (xo + ww, yo + crimp_h (image) + hn);
	    }
	}
    }

    /*
     * South East.
     */

    if (hs && we) {
	for (yo = 0; yo < hs; yo++) {
	    for (xo = 0; xo < we; xo++) {
		FILL_SE (xo + crimp_w (image) + ww, yo + crimp_h (image) + hn);
	    }
	}
    }

    /*
     * Central. Copy of the input image.
     */

    for (yo = hn, yi = 0; yi < crimp_h (image); yo++, yi++) {
	for (xo = ww, xi = 0; xi < crimp_w (image); xo++, xi++) {
	    COPY (xo, yo, xi, yi);
	}
    }

    Tcl_SetObjResult(interp, crimp_new_image_obj (result));
    return TCL_OK;

#!/bin/sh
# -*- tcl -*- \
exec tclsh "$0" ${1+"$@"}
# # ## ### ##### ######## ############# ######################
## Requisites

package require Tcl 8.5

# # ## ### ##### ######## ############# ######################
## Configuration

set simple_types {float grey32 grey16 grey8}
set multi_types  {rgb rgba fpcomplex}

array set accessor {
    float     {FLOATP v}
    grey32    {GREY32 v}
    grey16    {GREY16 v}
    grey8     {GREY8  v}
    hsv       {H h S s V v}
    rgb       {R r G g B b}
    rgba      {R r G g B b A a}
    fpcomplex {RE re IM im}
}

array set ctype {
    float     double
    grey32    {int   }
    grey16    {int   }
    grey8     {int   }
    hsv       {int   }
    rgb       {int   }
    rgba      {int   }
    fpcomplex double
}

# # ## ### ##### ######## ############# ######################
## Input

set basedir  [file dirname [file normalize [info script]]]
set template [apply {{template} {
    set c [open $template r]
    set data [read $c]
    close $c
    return $data
}} [file join $basedir binop_template.c]]

# # ## ### ##### ######## ############# ######################
## Generator core.

proc retrieve {type input image zero {expand 0} {alpha {}}} {
    global ctype accessor
    upvar 1 lines lines
    set variables {}

    set vprefix ${input}_
    set guard   in$input
    set ox      px$input
    set oy      py$input

    foreach {get var} $accessor($type) {
	lappend lines "$ctype($type) $vprefix$var = $guard ? $get ($image, $ox, $oy) : $zero;"
	lappend variables $vprefix$var
    }

    if {$expand} {
	set v [lindex $variables end]
	while {[llength $variables] < $expand} { lappend variables $v }
    }

    if {$alpha ne {}} {
	set var a
	lappend lines "$ctype($type) $vprefix$var = $guard ? OPAQUE : $zero;"
	lappend variables $vprefix$var
    }

    return $variables
}

proc assign {type avariables bvariables} {
    upvar 1 lines lines
    global accessor ctype

    lappend lines {}
    foreach av $avariables bv $bvariables {
	set zv z_[lindex [split $av _] end][lindex [split $bv _] end]
	lappend zvariables $zv
	lappend lines "$ctype($type) $zv = BINOP ($av, $bv);"
    }

    lappend lines {}
    foreach {put _} $accessor($type) zv $zvariables av $avariables bv $bvariables {
	lappend lines "$put (result, px, py) = BINOP_POST ($zv);"
    }
    return
}

proc assignA {type avariables bvariables} {
    upvar 1 lines lines
    global accessor

    # Match variable lists.
    while {[llength $bvariables] < [llength $avariables]} { lappend bvariables [lindex $bvariables end] }
    while {[llength $avariables] < [llength $bvariables]} { lappend avariables [lindex $avariables end] }

    lappend lines {}
    foreach {put _} $accessor($type) av $avariables bv $bvariables {
	if {$put eq "A"} {
	    lappend lines "$put (result, px, py) = MAX   ($av, $bv);"
	} else {
	    lappend lines "$put (result, px, py) = BINOP ($av, $bv);"
	}
    }
    return
}

proc gen {a b z} {
    global basedir template
    upvar 1 lines lines

    lappend map @TYPE_A@    $a
    lappend map @TYPE_B@    $b
    lappend map @TYPE_Z@    $z
    lappend map @TRANSFORM@ \t[join $lines \n\t]

    #puts \t[join $lines \n\t]

    set dst [file join $basedir binop_${a}_${b}_${z}.c]
    set ch  [open $dst w]
    puts -nonewline $ch [string map $map $template]
    close $ch
}

proc generate {a b z} {
    global accessor

    puts -nonewline "Generating $z = $a x $b ... "

    set la [dict size $accessor($a)]
    set lb [dict size $accessor($b)]
    set lz [dict size $accessor($z)]

    # .... generation with A channel ...

    if {($la == $lb) && ($la == $lz) &&
	[dict exists $accessor($a) A] &&
	[dict exists $accessor($b) A]} {
	puts {Av/A x Bv/A}
	# A, B vectors, identical length, element-wise operation.

	set av [retrieve $a a imageA BLACK]
	set bv [retrieve $b b imageB BLACK]
	assignA $z $av $bv
	gen $a $b $z
	return
    }

    if {($la > $lb) && ($la == $lz) &&
	[dict exists $accessor($a) A]} {
	puts {Av/A x B}
	# A vector, B scalar/vector expanded to match for element-wise
	# operation, pseudo-alpha for the scalar

	incr la -1
	set av [retrieve $a a imageA BLACK]
	set bv [retrieve $b b imageB BLACK $la alpha]
	assignA $z $av $bv
	gen $a $b $z
	return
    }

    if {($lb > $la) && ($lb == $lz) &&
	[dict exists $accessor($b) A]} {
	puts {A x Bv/A}
	# B vector, A scalar/vector expanded to match for element-wise
	# operation, pseudo-alpha for the scalar

	incr lb -1
	set av [retrieve $a a imageA BLACK $lb alpha]
	set bv [retrieve $b b imageB BLACK]
	assignA $z $av $bv
	gen $a $b $z
	return
    }

    # Generation without A channel.

    if {($la == $lb) && ($la == $lz)} {
	puts {Av x Bv}
	# A, B vectors, identical length, element-wise operation.

	set av [retrieve $a a imageA BLACK]
	set bv [retrieve $b b imageB BLACK]
	assign $z $av $bv
	gen $a $b $z
	return
    }

    if {($la > $lb) && ($la == $lz) && ($lb == 1)} {
	puts {Av x B}
	# A vector, B scalar expanded to match for element-wise operation

	set av [retrieve $a a imageA BLACK]
	set bv [retrieve $b b imageB BLACK $la]
	assign $z $av $bv
	gen $a $b $z
	return
    }

    if {($lb > $la) && ($lb == $lz) && ($la == 1)} {
	puts {A x Bv}
	# B vector, A scalar expanded to match for element-wise operation

	set av [retrieve $a a imageA BLACK $lb]
	set bv [retrieve $b b imageB BLACK]
	assign $z $av $bv
	gen $a $b $z
	return
    }

    puts "ERR"
    return -code error "Bad configuration ($z = $a x $b)"
}

# # ## ### ##### ######## ############# ######################
## Generate the various configurations

# # ## ### ##### ######## ############# ######################
## Simple vs simple.

generate float     float     float
generate float     grey16    float
generate float     grey32    float
generate float     grey8     float

generate grey16    float     float
generate grey16    grey16    grey16
generate grey16    grey32    grey32
generate grey16    grey8     grey16

generate grey32    float     float
generate grey32    grey16    grey32
generate grey32    grey32    grey32
generate grey32    grey8     grey32

generate grey8     float     float
generate grey8     grey16    grey16
generate grey8     grey32    grey32
generate grey8     grey8     grey8

# # ## ### ##### ######## ############# ######################
## Some operation always generate float, regardless of input types
## Examples: hypot, atan2, pow.

generate grey16    grey16    float
generate grey16    grey32    float
generate grey16    grey8     float

generate grey32    grey16    float
generate grey32    grey32    float
generate grey32    grey8     float

generate grey8     grey16    float
generate grey8     grey32    float
generate grey8     grey8     float

# # ## ### ##### ######## ############# ######################
## Thresholding against spatial floating point.

generate grey8     float     grey8
generate grey16    float     grey16
generate grey32    float     grey32

# # ## ### ##### ######## ############# ######################
## Complex vs. others, self and simple.

generate fpcomplex fpcomplex fpcomplex

generate fpcomplex float     fpcomplex
generate fpcomplex grey16    fpcomplex
generate fpcomplex grey32    fpcomplex
generate fpcomplex grey8     fpcomplex

generate float     fpcomplex fpcomplex
generate grey16    fpcomplex fpcomplex
generate grey32    fpcomplex fpcomplex
generate grey8     fpcomplex fpcomplex

# # ## ### ##### ######## ############# ######################
## RGB vs others, self and simple (grey8 only).

generate rgb       rgb       rgb

generate grey8     rgb       rgb
generate rgb       grey8     rgb

# # ## ### ##### ######## ############# ######################

generate hsv       hsv       hsv

# # ## ### ##### ######## ############# ######################
## RGBA vs self, RGB, and simple (grey8 only).

generate rgba  rgba  rgba

generate grey8 rgba  rgba
generate rgba  grey8 rgba

generate rgb   rgba  rgba
generate rgba  rgb   rgba

# # ## ### ##### ######## ############# ######################
exit

#!/bin/sh
# -*- tcl -*- \
exec tclsh "$0" ${1+"$@"}
# # ## ### ##### ######## ############# ######################
## Requisites

package require Tcl 8.5

# # ## ### ##### ######## ############# ######################
## Configuration

set simple_types {float grey32 grey16 grey8}
set multi_types  {rgb rgba fpcomplex}

array set accessor {
    float     {FLOATP v}
    grey32    {GREY32 v}
    grey16    {GREY16 v}
    grey8     {GREY8  v}
    hsv       {H h S s V v}
    rgb       {R r G g B b}
    rgba      {R r G g B b A a}
    fpcomplex {RE re IM im}
}

array set ctype {
    float     double
    grey32    int
    grey16    int
    grey8     int
    hsv       int
    rgb       int
    rgba      int
    fpcomplex double
}

# # ## ### ##### ######## ############# ######################
## Input

set basedir  [file dirname [file normalize [info script]]]
set template [apply {{template} {
    set c [open $template r]
    set data [read $c]
    close $c
    return $data
}} [file join $basedir unop_template.c]]

# # ## ### ##### ######## ############# ######################
## Generator core.

proc retrieve {type image zero {expand 0} {alpha {}}} {
    global ctype accessor
    upvar 1 lines lines
    set variables {}

    set vprefix _
    set guard   inside
    set ox      pxi
    set oy      pyi

    foreach {get var} $accessor($type) {
	lappend lines "$ctype($type) $vprefix$var = $guard ? $get ($image, $ox, $oy) : $zero;"
	lappend variables $vprefix$var
    }

    if {$expand} {
	set v [lindex $variables end]
	while {[llength $variables] < $expand} { lappend variables $v }
    }

    if {$alpha ne {}} {
	set var a
	lappend lines "$ctype($type) $vprefix$var = $guard ? OPAQUE : $zero;"
	lappend variables $vprefix$var
    }

    return $variables
}

proc assign {type avariables} {
    upvar 1 lines lines
    global accessor

    lappend lines {}
    foreach {put _} $accessor($type) av $avariables {
	lappend lines "$put (result, px, py) = UNOP ($av);"
    }
    return
}

proc assignA {type avariables} {
    upvar 1 lines lines
    global accessor

    lappend lines {}
    foreach {put _} $accessor($type) av $avariables {
	if {$put eq "A"} {
	    lappend lines "$put (result, px, py) = $av;"
	} else {
	    lappend lines "$put (result, px, py) = UNOP ($av);"
	}
    }
    return
}

proc gen {a z} {
    global basedir template
    upvar 1 lines lines

    lappend map @TYPE_IN@    $a
    lappend map @TYPE_OUT@   $z
    lappend map @TRANSFORM@ \t[join $lines \n\t]

    #puts \t[join $lines \n\t]

    set dst [file join $basedir unop_${a}.c]
    set ch  [open $dst w]
    puts -nonewline $ch [string map $map $template]
    close $ch
}

proc generate {a z} {
    global accessor

    puts -nonewline "Generating $z = $a ... "

    set la [dict size $accessor($a)]
    set lz [dict size $accessor($z)]

    # .... generation with A channel ...

    if {($la == $lz) &&
	[dict exists $accessor($a) A]} {
	puts {Av/A}
	# A vector, element-wise operation.

	set av [retrieve $a image BLACK]
	assignA $z $av
	gen $a $z
	return
    }

    # Generation without A channel.

    if {($la == $lz)} {
	puts {Av}
	# A vector, element-wise operation.

	set av [retrieve $a image BLACK]
	assign $z $av
	gen $a $z
	return
    }

    puts "ERR"
    return -code error "Bad configuration ($z = $a x $b)"
}

# # ## ### ##### ######## ############# ######################
## Generate the various configurations

# # ## ### ##### ######## ############# ######################

generate float     float
generate grey16    grey16
generate grey32    grey32
generate grey8     grey8

generate fpcomplex fpcomplex
generate rgb       rgb
generate rgba      rgba
generate hsv       hsv

# # ## ### ##### ######## ############# ######################
exit

 * #define OTYPE     unsigned char
 * (or whatever types are appropriate to the current instance)
 */

#define MAPNAME2(IN,OUT) crimp_piecewise_linear_map_ ## IN ## _ ## OUT
#define MAPNAME(IN,OUT) MAPNAME2(IN,OUT)

#define NEWNAME2(OUT) crimp_new_ ## OUT ## _at
#define NEWNAME(OUT) NEWNAME2(OUT)

#define STRINGIZE2(NAME) #NAME
#define STRINGIZE(NAME) STRINGIZE2(NAME)

crimp_image* inImage;
crimp_image* resultImage;

crimp_input(inImageObj, inImage, INTYPENAME);

resultImage = NEWNAME(OUTTYPENAME)(crimp_x (inImage), crimp_y (inImage),
				   crimp_w (inImage), crimp_h (inImage));

if (MAPNAME(INTYPENAME,OUTTYPENAME)(interp,
				    mapObj,
				    (size_t) (crimp_w (inImage) * crimp_h (inImage)),
		                    (const INTYPE*)inImage->pixel, 
                                    (size_t) 1,
				    (OUTTYPE*)resultImage->pixel, 
                                    (size_t) 1)
    != TCL_OK) {
    return TCL_ERROR;
}

/* x, y, w, h - Parameters of the output image. Provided by caller */

crimp_image*     result;
crimp_image*     image;

int px, py, lx, ly, ox, oy, pxi, pyi;

crimp_input (imageObj, image, float);

/*
 * Get the area of the input image to process.
 */

result = crimp_new_float_at (x, y, w, h);
ox = crimp_x (image);
oy = crimp_y (image);

/*
 * 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.
 */

for (py = 0, ly = y, pyi = y - oy;
     py < h;
     py++, ly++, pyi++) {

    for (px = 0, lx = x, pxi = x - ox;
	 px < w;
	 px++, lx++, pxi++) {

	int inside = crimp_inside (image, lx, ly);

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

	double _v = inside ? FLOATP (image, pxi, pyi) : BLACK;
	
	FLOATP (result, px, py) = UNOP (_v);
    }
}

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:
 */

/* x, y, w, h - Parameters of the output image. Provided by caller */

crimp_image*     result;
crimp_image*     image;

int px, py, lx, ly, ox, oy, pxi, pyi;

crimp_input (imageObj, image, fpcomplex);

/*
 * Get the area of the input image to process.
 */

result = crimp_new_fpcomplex_at (x, y, w, h);
ox = crimp_x (image);
oy = crimp_y (image);

/*
 * 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.
 */

for (py = 0, ly = y, pyi = y - oy;
     py < h;
     py++, ly++, pyi++) {

    for (px = 0, lx = x, pxi = x - ox;
	 px < w;
	 px++, lx++, pxi++) {

	int inside = crimp_inside (image, lx, ly);

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

	double _re = inside ? RE (image, pxi, pyi) : BLACK;
	double _im = inside ? IM (image, pxi, pyi) : BLACK;
	
	RE (result, px, py) = UNOP (_re);
	IM (result, px, py) = UNOP (_im);
    }
}

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:
 */

/* x, y, w, h - Parameters of the output image. Provided by caller */

crimp_image*     result;
crimp_image*     image;

int px, py, lx, ly, ox, oy, pxi, pyi;

crimp_input (imageObj, image, grey16);

/*
 * Get the area of the input image to process.
 */

result = crimp_new_grey16_at (x, y, w, h);
ox = crimp_x (image);
oy = crimp_y (image);

/*
 * 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.
 */

for (py = 0, ly = y, pyi = y - oy;
     py < h;
     py++, ly++, pyi++) {

    for (px = 0, lx = x, pxi = x - ox;
	 px < w;
	 px++, lx++, pxi++) {

	int inside = crimp_inside (image, lx, ly);

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

	int _v = inside ? GREY16 (image, pxi, pyi) : BLACK;
	
	GREY16 (result, px, py) = UNOP (_v);
    }
}

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:
 */

/* x, y, w, h - Parameters of the output image. Provided by caller */

crimp_image*     result;
crimp_image*     image;

int px, py, lx, ly, ox, oy, pxi, pyi;

crimp_input (imageObj, image, grey32);

/*
 * Get the area of the input image to process.
 */

result = crimp_new_grey32_at (x, y, w, h);
ox = crimp_x (image);
oy = crimp_y (image);

/*
 * 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.
 */

for (py = 0, ly = y, pyi = y - oy;
     py < h;
     py++, ly++, pyi++) {

    for (px = 0, lx = x, pxi = x - ox;
	 px < w;
	 px++, lx++, pxi++) {

	int inside = crimp_inside (image, lx, ly);

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

	int _v = inside ? GREY32 (image, pxi, pyi) : BLACK;
	
	GREY32 (result, px, py) = UNOP (_v);
    }
}

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:
 */

/* x, y, w, h - Parameters of the output image. Provided by caller */

crimp_image*     result;
crimp_image*     image;

int px, py, lx, ly, ox, oy, pxi, pyi;

crimp_input (imageObj, image, grey8);

/*
 * Get the area of the input image to process.
 */

result = crimp_new_grey8_at (x, y, w, h);
ox = crimp_x (image);
oy = crimp_y (image);

/*
 * 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.
 */

for (py = 0, ly = y, pyi = y - oy;
     py < h;
     py++, ly++, pyi++) {

    for (px = 0, lx = x, pxi = x - ox;
	 px < w;
	 px++, lx++, pxi++) {

	int inside = crimp_inside (image, lx, ly);

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

	int _v = inside ? GREY8 (image, pxi, pyi) : BLACK;
	
	GREY8 (result, px, py) = UNOP (_v);
    }
}

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:
 */