flot8.c

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

		p1[1] = p3[0] - p3[4];
		p1[3] = p3[1] - p3[5];
		p1[5] = p3[2] - p3[6];
		p1[7] = p3[3] - p3[7];
		p1 += 8;
		p3 += 8;
	}
    p1[0] = p3[0] * TWO;
    p1[2] = p3[1] * TWO;
    p1[4] = p3[2] * TWO;
    p1[6] = p3[3] * TWO;
    p1[1] = p1[3] = p1[5] = p1[7] = ZERO;

    /* Stages of IDCTs */

	p1 = y;
	p2 = y + xlen + 8;
    while (p1 < p2) {
        /* for LBT, scale first odd coefficient */
        if (type == LBT)
            p1[1] *= SQHALF;
        dctIII8(p1);
		p1 += 8;
	}

	/* Copy y back to x */
    q1 = x;
    p1 = y + 4;
    p3 = p1 + xlen;
    while (p1 < p3) {
        *q1++ = round(*p1++);
    }
}

/* Direct HLBT on a vector of blocks */

#define C00     0.9557930147983
#define S00     0.2940403252323
#define SQ1p8   7.07106781186548e-01 * 0.5

void hlbtvec(
        int16   *x,     /* vector of ints containing the i/o data */
        int     xlen,   /* length of vector (multiple of eight) */
        float   *y,     /* working vector, length = xlen + 4 */
        float   *z)     /* working vector, length = xlen + 4 */
{
	float   *p1, *p2, *p3, *p4;
    int16   *q1, *q2;
	int     i;

	/* Copy x onto y and fold borders */
    q1 = x;
    p1 = y + 2;
    p3 = p1 + xlen;
    while (p1 < p3) {
        *p1++ = *q1++;
    }
    p1 = p2 = y + 2;
	p3 = p4 = p1 + xlen;
    for (i = 2; i ; i--) {
		*(--p2) = *(p1++);
		*(p4++) = *(--p3);
	}

	/* Stages of DCTs */
	p1 = y;
    p2 = y + xlen + 4;
    while (p1 < p2) {
        dctII4(p1);
        p1[1] *= SQTWO;
        p1 += 4;
	}

	/* 1st stage of +1/-1 butterflies */
    z[0] = y[0];
    z[1] = y[2];
    p1 = z + 2;
    p2 = z + xlen - 2;
    p3 = y + 4;
    while(p1 < p2) {
		p1[0] = p3[0] + p3[1];
		p1[1] = p3[2] + p3[3];
        p1[2] = p3[0] - p3[1];
        p1[3] = p3[2] - p3[3];
        p1 += 4;
        p3 += 4;
	}
	p1[0] = p3[0];
	p1[1] = p3[2];

	/* 2nd stage of +1/-1 butterflies */
	p1 = y;
	p2 = y + xlen;
    p3 = z;
    while(p1 < p2) {
        p1[0] = p3[0] + p3[2];
        p1[1] = p3[1] + p3[3];
        p1[2] = p3[0] - p3[2];
        p1[3] = p3[1] - p3[3];
        p1 += 4;
        p3 += 4;
	}

	/* 3rd stage : plane rotations */
    p1 = z;
    p2 = z + xlen;
	p3 = y;
    while(p1 < p2) {
        p1[0] = p3[0];
        p1[2] = p3[1];
        p1[1] = p3[2] * C00 - p3[3] * S00;
        p1[3] = p3[2] * S00 + p3[3] * C00;
        p1 += 4;
        p3 += 4;
	}

    /* 4th stage : combines DCs of adjacent blocks */
    q1 = x;
    q2 = x + xlen;
    p3 = z;
    while(q1 < q2) {
        q1[0] = round(SQ1p8 * (p3[0] + p3[4]));
        q1[1] = round(SQ1p8 * (p3[0] - p3[4]));
        q1[2] = round(HALF * p3[1]);
        q1[4] = round(HALF * p3[2]);
        q1[6] = round(HALF * p3[3]);
        q1[3] = round(HALF * p3[5]);
        q1[5] = round(HALF * p3[6]);
        q1[7] = round(HALF * p3[7]);
        q1 += 8;
        p3 += 8;
    }
}

/* Inverse HLBT on a vector of blocks */

void ihlbtvec(
        int16   *x,     /* vector of ints containing the i/o data */
        int     xlen,   /* length of vector (multiple of eight) */
        float   *y,     /* working vector, length = xlen + 4 */
        float   *z)     /* working vector, length = xlen + 4 */
{
	float   *p1, *p2, *p3;
    int16   *q1, *q3;

    /* 4th stage : combines DCs of adjacent blocks */
    p1 = z;
    p2 = z + xlen;
    q3 = x;
    while(p1 < p2) {
        p1[0] = SQ1p8 * (q3[0] + (float) q3[1]);
        p1[4] = SQ1p8 * (q3[0] - (float) q3[1]);
        p1[1] = HALF * q3[2];
        p1[2] = HALF * q3[4];
        p1[3] = HALF * q3[6];
        p1[5] = HALF * q3[3];
        p1[6] = HALF * q3[5];
        p1[7] = HALF * q3[7];
        p1 += 8;
        q3 += 8;
    }

    /* 3rd stage : plane rotations */
	p1 = y;
	p2 = y + xlen;
    p3 = z;
    while(p1 < p2) {
        p1[0] = p3[0];
        p1[1] = p3[2];
        p1[2] = p3[1] * C00 + p3[3] * S00;
        p1[3] = p3[3] * C00 - p3[1] * S00;
        p1 += 4;
        p3 += 4;
	}

	/* 2nd stage of +1/-1 butterflies */
    p1 = z;
    p2 = z + xlen;
	p3 = y;
    while(p1 < p2) {
        p1[0] = p3[0] + p3[2];
        p1[1] = p3[1] + p3[3];
        p1[2] = p3[0] - p3[2];
        p1[3] = p3[1] - p3[3];
        p1 += 4;
        p3 += 4;
	}

	/* 1st stage of +1/-1 butterflies */
    y[0] = z[0] * TWO;
    y[2] = z[1] * TWO;
    y[1] = y[3] = ZERO;
    p1 = y + 4;
	p2 = y + xlen;
    p3 = z + 2;
    while(p1 < p2) {
        p1[0] = p3[0] + p3[2];
        p1[2] = p3[1] + p3[3];
        p1[1] = p3[0] - p3[2];
        p1[3] = p3[1] - p3[3];
        p1 += 4;
        p3 += 4;
	}
    p1[0] = p3[0] * TWO;
    p1[2] = p3[1] * TWO;
    p1[1] = p1[3] = ZERO;

    /* Stages of IDCTs */
	p1 = y;
    p2 = y + xlen + 4;
    while (p1 < p2) {
        p1[1] *= SQHALF;
        dctIII4(p1);
        p1 += 4;
	}

	/* Copy y back to x */
    q1 = x;
    p1 = y + 2;
    p3 = p1 + xlen;
    while (p1 < p3) {
        *q1++ = round(*p1++);
    }
}

/* --------- PART 2:  IMAGE TRANSFORM ROUTINES --------- */

/* Compute the forward DCT of an image */

void forward_dct(pixel **p, int16 **img, int nrows, int ncols)
{
	int     k, j;
	int16   *t;

    /* Allocate memory for auxiliary vector */
    t = allocint(nrows);

    /* Transform rows */
    for (k = 0 ; k < nrows ; k++) {
        for (j = 0 ; j < ncols; j++) {
			/* remove mean and scale up by a factor of 8 */
            img[k][j] = (((int16) p[k][j]) - 128) << 3;
        }
        dctvec(img[k], ncols);
	}

	/* Transform columns */
    for (k = 0 ; k < ncols ; k++) {
        for (j = 0 ; j < nrows ; j++)
			t[j] = img[j][k];
        dctvec(t, nrows);
        for (j = 0 ; j < nrows ; j++)
			img[j][k] = t[j];
	}

    /* Free auxiliary vector */
	free(t);
}

/* Compute the inverse DCT of an image */

void inverse_dct(int16 **img, pixel **p, int nrows, int ncols)
{
	int     k, j;
    int16   *t, v;

    /* Allocate memory for auxiliary vector */
    t = allocint(nrows);

	/* Transform columns */
    for (k = 0 ; k < ncols ; k++) {
        for (j = 0 ; j < nrows ; j++)
			t[j] = img[j][k];
        idctvec(t, nrows);
        for (j = 0 ; j < nrows ; j++)
			img[j][k] = t[j];
	}

    /* Transform rows */
    for (k = 0 ; k < nrows ; k++) {
		idctvec(img[k], ncols);
        for (j = 0 ; j < ncols; j++) {
            v = img[k][j];
            if (v < -1024) v = -1024;
            if (v > 1023) v = 1023;
			/* scale down by a factor of 8 & add back mean */
            p[k][j] = (v >> 3) + 128;
        }
	}

    /* Free auxiliary vector */
	free(t);
}


/* Compute the forward LOT of an image */

void forward_lot(pixel **p, int16 **img, int nrows, int ncols, int type)
{
    int     k, j, nf;
	int16   *t;
    float   *y, *z;

    /* Allocate memory for auxiliary vectors */
    t = allocint(nrows);
    nf = nrows;
    if (ncols > nf) nf = ncols;
    nf += 8;
    y = allocfloat(nf);
    z = allocfloat(nf);

    /* Transform rows */
    for (k = 0 ; k < nrows ; k++) {
        for (j = 0 ; j < ncols; j++) {
			/* remove mean and scale up by a factor of 8 */
            img[k][j] = (((int16) p[k][j]) - 128) << 3;
        }
        lotvec(img[k], ncols, y, z, type);
	}

	/* Transform columns */
    for (k = 0 ; k < ncols ; k++) {
        for (j = 0 ; j < nrows ; j++)
			t[j] = img[j][k];
        lotvec(t, nrows, y, z, type);
        for (j = 0 ; j < nrows ; j++)
			img[j][k] = t[j];
	}

    /* Free auxiliary vector */
	free(t);
    free(y);
    free(z);
}

/* Compute the inverse LOT of an image */

void inverse_lot(int16 **img, pixel **p, int nrows, int ncols, int type)
{
    int     k, j, nf;
    int16   *t, v;
    float   *y, *z;

    /* Allocate memory for auxiliary vector */
    t = allocint(nrows);
    nf = nrows;
    if (ncols > nf) nf = ncols;
    nf += 8;
    y = allocfloat(nf);
    z = allocfloat(nf);

	/* Transform columns */
    for (k = 0 ; k < ncols ; k++) {
        for (j = 0 ; j < nrows ; j++)
			t[j] = img[j][k];
        ilotvec(t, nrows, y, z, type);
        for (j = 0 ; j < nrows ; j++)
			img[j][k] = t[j];
	}

    /* Transform rows */
    for (k = 0 ; k < nrows ; k++) {
        ilotvec(img[k], ncols, y, z, type);
        for (j = 0 ; j < ncols; j++) {
            v = img[k][j];
            if (v < -1024) v = -1024;
            if (v > 1023) v = 1023;
			/* scale down by a factor of 8 & add back mean */
            p[k][j] = (v >> 3) + 128;
        }
	}

    /* Free auxiliary vector */
	free(t);
}

/* Compute the forward HLBT of an image */

void forward_hlbt(pixel **p, int16 **img, int nrows, int ncols)
{
    int     k, j, nf;
	int16   *t;
    float   *y, *z;

    /* Allocate memory for auxiliary vectors */
    t = allocint(nrows);
    nf = nrows;
    if (ncols > nf) nf = ncols;
    nf += 4;
    y = allocfloat(nf);
    z = allocfloat(nf);

    /* Transform rows */
    for (k = 0 ; k < nrows ; k++) {
        for (j = 0 ; j < ncols; j++) {
			/* remove mean and scale up by a factor of 8 */
            img[k][j] = (((int16) p[k][j]) - 128) << 3;
        }
        hlbtvec(img[k], ncols, y, z);
	}

	/* Transform columns */
    for (k = 0 ; k < ncols ; k++) {
        for (j = 0 ; j < nrows ; j++)
			t[j] = img[j][k];
        hlbtvec(t, nrows, y, z);
        for (j = 0 ; j < nrows ; j++)
			img[j][k] = t[j];
	}

    /* Free auxiliary vector */
	free(t);
    free(y);
    free(z);
}

/* Compute the inverse HLBT of an image */

void inverse_hlbt(int16 **img, pixel **p, int nrows, int ncols)
{
    int     k, j, nf;
    int16   *t, v;
    float   *y, *z;

    /* Allocate memory for auxiliary vector */
    t = allocint(nrows);
    nf = nrows;
    if (ncols > nf) nf = ncols;
    nf += 8;
    y = allocfloat(nf);
    z = allocfloat(nf);

	/* Transform columns */
    for (k = 0 ; k < ncols ; k++) {
        for (j = 0 ; j < nrows ; j++)
			t[j] = img[j][k];
        ihlbtvec(t, nrows, y, z);
        for (j = 0 ; j < nrows ; j++)
			img[j][k] = t[j];
	}

    /* Transform rows */
    for (k = 0 ; k < nrows ; k++) {
        ihlbtvec(img[k], ncols, y, z);
        for (j = 0 ; j < ncols; j++) {
            v = img[k][j];
            if (v < -1024) v = -1024;
            if (v > 1023) v = 1023;
			/* scale down by a factor of 8 & add back mean */
            p[k][j] = (v >> 3) + 128;
        }
	}

    /* Free auxiliary vector */
	free(t);
}