}
	    }
	}
    }

    return [map_$type $image {*}$args]
}

proc ::crimp::downsample {image factor} {
    set type [TypeOf $image]
    set f downsample_$type
    if {![Has $f]} {
	return -code error "Unable to downsample images of type \"$type\""
    }

    return [$f $image $factor]
}

proc ::crimp::upsample {image factor} {
    set type [TypeOf $image]
    set f upsample_$type
    if {![Has $f]} {
	return -code error "Unable to upsample images of type \"$type\""
    }

    return [$f $image $factor]
}

proc ::crimp::decimate {image factor kernel} {
    # Combines downsampling with a pre-processing step applying a
    # low-pass filter to avoid aliasing of higher image frequencies.

    # We assume that the low-pass filter is separable, and the kernel
    # is the 1-D horizontal form of it. We compute the vertical form
    # on our own, transposing the kernel.

    # NOTE: This implementation, while easy conceptually, is not very
    # efficient, because it does the filtering on the input image,
    # before downsampling.

    # FUTURE: Write C level primitive integrating filter and sampler,
    # computing the filter only for the pixels which go into the
    # result.

    return [downsample \
		[convolve $image $kernel [kernel transpose $kernel]] \
		$factor]
}

proc ::crimp::interpolate {image factor kernel} {
    # Combines upsampling with a post-processing step applying a
    # low-pass filter to copies of the image at higher image
    # frequencies.

    # We assume that the low-pass filter is separable, and the kernel
    # is the 1-D horizontal form of it. We compute the vertical form
    # on our own, transposing the kernel.

    # NOTE: This implementation, while easy conceptually, is not very
    # efficient, because it does the filtering on the full output image,
    # after upsampling.

    # FUTURE: Write C level primitive integrating filter and sampler,
    # computing the filter only for the actually new pixels, and use
    # polyphase restructuring.

    # DANGER: This assumes that the filter, applied to the original
    # pixels leaves them untouched. I.e. scaled center weight is 1.
    # The easy implementation here does not have this assumption.

    return [convolve [upsample $image $factor] \
		$kernel [kernel transpose $kernel]]
}

proc ::crimp::split {image} {
    set type [TypeOf $image]
    if {![Has split_$type]} {
	return -code error "Unable to split images of type \"$type\""
    }
    return [split_$type $image]

namespace eval ::crimp {
    namespace export type width height dimensions channels
    namespace export read write convert join flip split table
    namespace export invert solarize gamma degamma remap map
    namespace export wavy psychedelia matrix blend over blank
    namespace export setalpha histogram max min screen add
    namespace export subtract difference multiply convolve
    namespace export downsample upsample decimate interpolate
    namespace export kernel expand
    #
    namespace ensemble create
}

# # ## ### ##### ######## #############
return

def op_convolve_pseudoedges {
    label {Pseudo Edges}
    setup {
	variable K [crimp kernel make {
	    {1  2  4  2 1}
	    {2  4  8  4 2}
	    {4  8 16  8 4}
	    {2  4  8  4 2}
	    {1  2  4  2 1}}]

	# Separable kernel, compute the horizontal and vertical kernels.
	variable Kx [crimp kernel make {{1 2 4 2 1}}]
	variable Ky [crimp kernel transpose $Kx]
    }
    setup_image {
	# show_image [crimp convolve [base] $K]
	# Separable kernel, convolve x and y separately. Same result
	# as for the combined kernel, but faster.
	set n [crimp add [base] [crimp difference [base] [crimp convolve [base] $Kx $Ky]]]
	set m [crimp blank grey8 {*}[crimp dimensions $n] 255]
	show_image [crimp setalpha $n $m]
    }
}

def op_decimate2 {
    label Decimate\u21932
    setup {
	set K [crimp kernel make {{1 2 1}}]
    }
    setup_image {
	set n [crimp decimate [base] 2 $K]
	set m [crimp blank grey8 {*}[crimp dimensions $n] 255]
	show_image [crimp setalpha $n $m]
    }
}

def op_decimate4 {
    label Decimate\u21934
    setup {
	set K [crimp kernel make {{1 2 1}}]
    }
    setup_image {
	set n [crimp decimate \
		   [crimp decimate [base] 2 $K] \
		   2 $K]
	set m [crimp blank grey8 {*}[crimp dimensions $n] 255]
	show_image [crimp setalpha $n $m]
    }
}

def op_decint2 {
    label "Decimate\u21932, Interpolate\u21912"
    setup {
	set KD [crimp kernel make {{1 2 1}}]
	set KI [crimp kernel make {{1 2 1}} 2]
    }
    setup_image {
	set n [crimp interpolate \
		   [crimp decimate [base] 2 $KD] \
		   2 $KI]
	set m [crimp blank grey8 {*}[crimp dimensions $n] 255]
	show_image [crimp setalpha $n $m]
    }
}

def op_downsample2 {
    label Downsample\u21932
    setup_image {
	show_image [crimp downsample [base] 2]
    }
}

def op_downsample3 {
    label Downsample\u21933
    setup_image {
	show_image [crimp downsample [base] 3]
    }
}

def op_downsample4 {
    label Downsample\u21934
    setup_image {
	show_image [crimp::downsample_rgba [base] 4]
    }
}

def op_interpolate2 {
    label Interpolate\u21912
    setup {
	# Tent kernel.
	set K [crimp kernel make {{1 2 1}} 2]
    }
    setup_image {
	set n [crimp interpolate [base] 2 $K]
	set m [crimp blank grey8 {*}[crimp dimensions $n] 255]
	show_image [crimp setalpha $n $m]
    }
}

def op_interpolate4 {
    label Interpolate\u21914
    setup {
	# Tent kernel.
	set K [crimp kernel make {{1 2 1}} 2]
    }
    setup_image {
	set n [crimp interpolate \
		   [crimp interpolate [base] 2 $K] \
		   2 $K]
	set m [crimp blank grey8 {*}[crimp dimensions $n] 255]
	show_image [crimp setalpha $n $m]
    }
}

def op_upsample2 {
    label Upsample\u21912
    setup_image {
	show_image [crimp upsample [base] 2]
    }
}

def op_upsample3 {
    label Upsample\u21913
    setup_image {
	show_image [crimp upsample [base] 3]
    }
}

def op_upsample4 {
    label Upsample\u21914
    setup_image {
	show_image [crimp upsample [base] 4]
    }
}

downsample_grey8
Tcl_Obj* imageObj
int      factor

/*
 * The input image is downsampled by storing only every 'factor' pixel into
 * the result. Note that this method of shrinking an image causes image
 * frequencies above the nyquist threshold of the result to be aliased into
 * the range.
 *
 * The input image has to be convolved with a low-pass filter first, to avoid
 * such artefacts. The integrated combination of such a filter with
 * downsampling is called 'decimation'. This is but one step in the generation
 * of image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi;

crimp_input (imageObj, image, grey8);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w/factor, image->h/factor);

for (yo = 0, yi = 0; yo < result->h; yo++, yi += factor) {
    for (xo = 0, xi = 0; xo < result->w; xo++, xi += factor) {

	GREY8 (result, xo, yo) = GREY8 (image, xi, yi);
    }
}

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

downsample_hsv
Tcl_Obj* imageObj
int      factor

/*
 * The input image is downsampled by storing only every 'factor' pixel into
 * the result. Note that this method of shrinking an image causes image
 * frequencies above the nyquist threshold of the result to be aliased into
 * the range.
 *
 * The input image has to be convolved with a low-pass filter first, to avoid
 * such artefacts. The integrated combination of such a filter with
 * downsampling is called 'decimation'. This is but one step in the generation
 * of image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi;

crimp_input (imageObj, image, hsv);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w/factor, image->h/factor);

for (yo = 0, yi = 0; yo < result->h; yo++, yi += factor) {
    for (xo = 0, xi = 0; xo < result->w; xo++, xi += factor) {

	H (result, xo, yo) = H (image, xi, yi);
	S (result, xo, yo) = S (image, xi, yi);
	V (result, xo, yo) = V (image, xi, yi);
    }
}

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

downsample_rgb
Tcl_Obj* imageObj
int      factor

/*
 * The input image is downsampled by storing only every 'factor' pixel into
 * the result. Note that this method of shrinking an image causes image
 * frequencies above the nyquist threshold of the result to be aliased into
 * the range.
 *
 * The input image has to be convolved with a low-pass filter first, to avoid
 * such artefacts. The integrated combination of such a filter with
 * downsampling is called 'decimation'. This is but one step in the generation
 * of image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi;

crimp_input (imageObj, image, rgb);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w/factor, image->h/factor);

for (yo = 0, yi = 0; yo < result->h; yo++, yi += factor) {
    for (xo = 0, xi = 0; xo < result->w; xo++, xi += factor) {

	R (result, xo, yo) = R (image, xi, yi);
	G (result, xo, yo) = G (image, xi, yi);
	B (result, xo, yo) = B (image, xi, yi);
    }
}

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

downsample_rgba
Tcl_Obj* imageObj
int      factor

/*
 * The input image is downsampled by storing only every 'factor' pixel into
 * the result. Note that this method of shrinking an image causes image
 * frequencies above the nyquist threshold of the result to be aliased into
 * the range.
 *
 * The input image has to be convolved with a low-pass filter first, to avoid
 * such artefacts. The integrated combination of such a filter with
 * downsampling is called 'decimation'. This is but one step in the generation
 * of image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi;

crimp_input (imageObj, image, rgba);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w/factor, image->h/factor);

for (yo = 0, yi = 0; yo < result->h; yo++, yi += factor) {
    for (xo = 0, xi = 0; xo < result->w; xo++, xi += factor) {

	R (result, xo, yo) = R (image, xi, yi);
	G (result, xo, yo) = G (image, xi, yi);
	B (result, xo, yo) = B (image, xi, yi);
	A (result, xo, yo) = A (image, xi, yi);
    }
}

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

upsample_grey8
Tcl_Obj* imageObj
int      factor

/*
 * The input image is upsampled by inserting 'factor-1' 0-pixels after every
 * pixel of the input. Note that this method of expanding an image introduces
 * copies of the input to appear at higher frequencies.
 *
 * The output image has to be convolved with a low-pass filter after expansion
 * to avoid such artefacts. The integrated combination of upsampling and such
 * a filter is called 'interpolation'. This is but one step in the generation
 * of difference image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi, dx, dy;

crimp_input (imageObj, image, grey8);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w*factor, image->h*factor);

for (yo = 0, yi = 0; yi < image->h; yo += factor, yi ++) {
    for (xo = 0, xi = 0; xi < image->w; xo += factor, xi ++) {

	/* Copy the pixel */
	GREY8 (result, xo, yo) = GREY8 (image, xi, yi);

	/* And insert factor black (0) pixels after */
	for (dx = 1; dx < factor; dx++) {
	    GREY8 (result, xo + dx, yo) = BLACK;
	}
    }

    /* And insert factor black lines after the intput line*/
    for (dy = 1; dy < factor; dy++) {
	for (xo = 0; xo < result->w; xo++) {
	    GREY8 (result, xo, yo + dy) = BLACK;
	}
    }
}

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

upsample_hsv
Tcl_Obj* imageObj
int      factor

/*
 * The input image is upsampled by inserting 'factor-1' 0-pixels after every
 * pixel of the input. Note that this method of expanding an image introduces
 * copies of the input to appear at higher frequencies.
 *
 * The output image has to be convolved with a low-pass filter after expansion
 * to avoid such artefacts. The integrated combination of upsampling and such
 * a filter is called 'interpolation'. This is but one step in the generation
 * of difference image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi, dx, dy;

crimp_input (imageObj, image, hsv);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w*factor, image->h*factor);

for (yo = 0, yi = 0; yi < image->h; yo += factor, yi ++) {
    for (xo = 0, xi = 0; xi < image->w; xo += factor, xi ++) {

	/* Copy the pixel */
	H (result, xo, yo) = H (image, xi, yi);
	S (result, xo, yo) = S (image, xi, yi);
	V (result, xo, yo) = V (image, xi, yi);

	/* And insert factor black (0) pixels after */
	for (dx = 1; dx < factor; dx++) {
	    H (result, xo + dx, yo) = BLACK;
	    S (result, xo + dx, yo) = BLACK;
	    V (result, xo + dx, yo) = BLACK;
	}
    }

    /* And insert factor black lines after the intput line*/
    for (dy = 1; dy < factor; dy++) {
	for (xo = 0; xo < result->w; xo++) {
	    H (result, xo, yo + dy) = BLACK;
	    S (result, xo, yo + dy) = BLACK;
	    V (result, xo, yo + dy) = BLACK;
	}
    }
}

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

upsample_rgb
Tcl_Obj* imageObj
int      factor

/*
 * The input image is upsampled by inserting 'factor-1' 0-pixels after every
 * pixel of the input. Note that this method of expanding an image introduces
 * copies of the input to appear at higher frequencies.
 *
 * The output image has to be convolved with a low-pass filter after expansion
 * to avoid such artefacts. The integrated combination of upsampling and such
 * a filter is called 'interpolation'. This is but one step in the generation
 * of difference image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi, dx, dy;

crimp_input (imageObj, image, rgb);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w*factor, image->h*factor);

for (yo = 0, yi = 0; yi < image->h; yo += factor, yi ++) {
    for (xo = 0, xi = 0; xi < image->w; xo += factor, xi ++) {

	/* Copy the pixel */
	R (result, xo, yo) = R (image, xi, yi);
	G (result, xo, yo) = G (image, xi, yi);
	B (result, xo, yo) = B (image, xi, yi);

	/* And insert factor black (0) pixels after */
	for (dx = 1; dx < factor; dx++) {
	    R (result, xo + dx, yo) = BLACK;
	    G (result, xo + dx, yo) = BLACK;
	    B (result, xo + dx, yo) = BLACK;
	}
    }

    /* And insert factor black lines after the intput line*/
    for (dy = 1; dy < factor; dy++) {
	for (xo = 0; xo < result->w; xo++) {
	    R (result, xo, yo + dy) = BLACK;
	    G (result, xo, yo + dy) = BLACK;
	    B (result, xo, yo + dy) = BLACK;
	}
    }
}

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

upsample_rgba
Tcl_Obj* imageObj
int      factor

/*
 * The input image is upsampled by inserting 'factor-1' 0-pixels after every
 * pixel of the input. Note that this method of expanding an image introduces
 * copies of the input to appear at higher frequencies.
 *
 * The output image has to be convolved with a low-pass filter after expansion
 * to avoid such artefacts. The integrated combination of upsampling and such
 * a filter is called 'interpolation'. This is but one step in the generation
 * of difference image pyramids.
 */

crimp_image* image;
crimp_image* result;
int          xo, yo, xi, yi, dx, dy;

crimp_input (imageObj, image, rgba);
if (factor < 1) {
    Tcl_SetResult(interp, "bad sampling factor, expected integer > 0", TCL_STATIC);
    return TCL_ERROR;
}

if (factor == 1) {
    Tcl_SetObjResult(interp, imageObj);
    return TCL_OK;
}

result = crimp_new (image->itype, image->w*factor, image->h*factor);

for (yo = 0, yi = 0; yi < image->h; yo += factor, yi ++) {
    for (xo = 0, xi = 0; xi < image->w; xo += factor, xi ++) {

	/* Copy the pixel */
	R (result, xo, yo) = R (image, xi, yi);
	G (result, xo, yo) = G (image, xi, yi);
	B (result, xo, yo) = B (image, xi, yi);
	A (result, xo, yo) = A (image, xi, yi);

	/* And insert factor black (0) pixels after */
	for (dx = 1; dx < factor; dx++) {
	    R (result, xo + dx, yo) = BLACK;
	    G (result, xo + dx, yo) = BLACK;
	    B (result, xo + dx, yo) = BLACK;
	    A (result, xo + dx, yo) = OPAQUE;
	}
    }

    /* And insert factor black lines after the intput line*/
    for (dy = 1; dy < factor; dy++) {
	for (xo = 0; xo < result->w; xo++) {
	    R (result, xo, yo + dy) = BLACK;
	    G (result, xo, yo + dy) = BLACK;
	    B (result, xo, yo + dy) = BLACK;
	    A (result, xo, yo + dy) = OPAQUE;
	}
    }
}

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