testlot.c

lapped transforms for images. lapped transform is a very fast algorithm to replace dct.

/* 
 * Test program for routines for fast DCT and lapped transforms for images
 *
 * H. S. Malvar, June '97
 *
 * Refs: - Malvar, Signal Processing with Lapped Transforms, book, 1992.
 *       - Malvar, LBT/HLBT paper, IEEE Trans. SP, Apr. 1998.
 *
 * History:
 *
 * 06/97 - HSM, first version
 * 10/00 - HSM, public release
 *
 * Disclaimer: use these programs at your own risk; no guarantees!
 * The code is reasonably fast but not optimized for memory allocation or speed.
 *
 * Have fun with lapped transforms!
 *
 * PS1: this software is (c) 1993-2000 H. S. Malvar, all rights reserved.
 * You may use this code or derivatives at will, but please quote the source ...
 *
 * PS2: this program is "postcardware"; if you decide to use it I do charge a
 * license fee: you should send me a postcard of the city in which you work :-)
 * please send it to: H. Malvar, Microsoft Research, One Microsoft Way,
 * Redmond, WA 98052, USA.
 */

#include <stdlib.h>
#include <stdio.h>
#include "flot8.h"

/* Allocate a block of memory */

void * alloc(uint nitems, uint size)
{
	void	*retval;

	if (!(retval = calloc(nitems, size))) {
		printf("Not enough memory\n");
		exit(1);
	}
	return retval;
}

/* Read & write raw imge data */

void readimg(char *s, pixel **x, int nrows, int ncols)
{
	FILE	*f;
	int		i, j;

	/* try to open file for reading */
	if ((f = fopen(s, "rb")) == NULL) {
		printf("File %s not found\n", s);
		exit (1);
	}

	/* read the rows */
	for (i = 0; i < nrows; i++) {

		/* check if sizes match during reading */
		if (fread(x[i], sizeof (pixel), ncols, f) != ncols) {
			printf("Error reading file %s\n", s);
			exit (1);
		}
	}

	fclose(f);
}


void writeimg(char *s, pixel **x, int nrows, int ncols)
{
	FILE	*f;
	int		i, j, k, delta, n;

	/* try to open file */
	if ((f = fopen(s, "wb")) == NULL) {
		printf("Error creating file %s\n", s);
		exit(1);
	}

	/* write image */
	for (i = 0; i < nrows; i++) {

		/* write the rows */
		if (fwrite (x[i], sizeof (pixel), ncols, f) != ncols) {
			printf ("Error writing file %s (disk full?)\n", s);
			exit (1);
		}
	}

	fclose (f);
}


/* Test program: very crude ... :-)
   - read input image file
   - perform forward transform
   - quantize all transform coefficients
   - perform inverse transform
   - write output image file
*/

void help(void)
{
	printf(
		"\nSyntax: testlot input_image nrows ncols type qstep output_image\n"
		"\n"
		"Input and output images are in raw format, grayscale, 8 bits/pixel\n"
		"(e.g. PaintShop Pro raw format).\n"
		"\n"
		"  nrows and ncols must be a multiple of the blocksize = 8\n"
		"  type specifies the transform: 1 = DCT, 2 = LOT, 3 = LBT, 4 = HLBT\n"
		"  qstep specifies the quantization step size\n"
		"\n"
		"Example: testlot bikes.raw 512 768 1 150 new.raw\n"
		"\n");
	exit(0);
}

void main(int argc, char *argv[])
{
	pixel	**p;	/* pixel data */
	int16	**img;	/* transformed image data */
	int16	qv;
	int		i, j, nrows, ncols, type;
	float	step, invstep;

	if (argc < 7) help();

	/* define image size */
	nrows = atoi(argv[2]);
	ncols = atoi(argv[3]);

	/* define transform type */
	type = atoi(argv[4]);
	if ((type < DCT) || (type > HLBT)) {
		printf("\nOops: invalid transform type.\n");
		help();
	}

	/* define quantization step size */
	step = atoi(argv[5]);
	invstep = 1.0 / step;

	/* allocate memory for pixel and image data */
	p = alloc(nrows, sizeof(pixel *));
	img = alloc(nrows, sizeof(int16 *));
	for (i = 0; i < nrows; i++) {
		/* alloc mem for the row */
		p[i] = alloc(ncols, sizeof(pixel));
		img[i] = alloc(ncols, sizeof(int16));
	}

	/* read pixels */
	readimg(argv[1], p, nrows, ncols);

	/* compute forward DCT, LOT, LBT, or HLBT */
	switch (type) {
	case DCT:
		forward_dct(p, img, nrows, ncols);
		break;
	case LOT:
	case LBT:
		forward_lot(p, img, nrows, ncols, type);
		break;
	case HLBT:
		forward_hlbt(p, img, nrows, ncols);
	}

	/*  quantize transform coefficients.
		Dynamic range of coefficients is |img| <= 2048, since the forward
		transform routines have a gain of 8.
		Add code for some good entropy coder for the quantized values qv
		below, and you'll have a built a transform coder (see if you can beat
		JPEG - I bet you can!)
	*/
	for (i = 0; i < nrows; i++) {
		for (j = 0; j < ncols; j++) {
			if (img[i][j] >= 0) {
				qv = (int16) (img[i][j] * invstep + 0.5);
			} else {
				qv = (int16) (img[i][j] * invstep - 0.5);
			}
			img[i][j]  = step * qv;
		}
	}

	/* compute inverse DCT, LOT, LBT, or HLBT */
	switch (type) {
	case DCT:
		inverse_dct(img, p, nrows, ncols);
		break;
	case LOT:
	case LBT:
		inverse_lot(img, p, nrows, ncols, type);
		break;
	case HLBT:
		inverse_hlbt(img, p, nrows, ncols);
	}

	/* write pixels */
	writeimg(argv[6], p, nrows, ncols);

	/* free memory */
	for (i = 0; i < nrows; i++) {
		free(img[i]); free(p[i]);
	}
	free(img); free(p);
}