spt.c

用C++语言实现的基于小波分析的源代码,实现了小波分析的诸多算法

/****************

I unrolled some if's on 7-2 : did it slow this down? !?

-----

todos for speed :

	1. don't do the "S" and "P" in separate steps : inline them and
		thereby eliminate the need for a temp workspace as well
	2. don't actually transpose; include a "step" in the _tdec_ call

******************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <crblib/inc.h>

/** line is just a temp workspace **/

static void do_tdec_line(int *data,int *line,int len)
{
int x,a,b,E,W,ND,pred;
int *high;

//  if ( len <= 2 ) errexit("special-cased to len > 2 , scumbag");
//	if ( len&1 ) errexit("len shouldn't be odd");

	for(x=0;x<len;x += 2) {
		a = data[x];
		b = data[x+1];

		line[x] = (a+b)>>1;
		line[x+1] = a-b;
	}

	line[1] += (( line[2] - line[0] + line[3] ) >> 1);

	for(x=2;x<=(len-4);x += 2) {
		E = line[x + 2] - line[x - 2];
		W = line[x + 2] - line[x];
		ND = line[x + 3];
		pred = ( E + E + W + ND + ND ) >> 3;
		line[x+1] += pred;
	}

	// x = len-2;
	line[len-1] += line[len-2] - line[len-4];

	/** line is now interleaved LHLHLH ; reshuffle to LLLLHHHH **/

	high = data + (len>>1);
	for(x=0;x<(len>>1);x++) {
		data[x] = line[x+x];
		high[x] = line[x+x+1];
	}
}

static void un_tdec_line(int *data,int *line,int len)
{
int x,a,A,D;
int E,W,ND,pred;
int *high;

	/** data is in LLLHHHH order ; we must deshuffle it into line ***/

	high = data + (len>>1);
	for(x=0;x<(len>>1);x++) {
		line[x+x] = data[x];
		line[x+x+1] = high[x];
	}	

	//** run backwards cuz of ND

	// x = len-2;
	line[len-1] -= line[len-2] - line[len-4];

	for(x=len-4;x>=2;x -= 2) {
		E = line[x + 2] - line[x - 2];
		W = line[x + 2] - line[x];
		ND = line[x + 3];
		pred = ( E + E + W + ND + ND ) >> 3;
		line[x+1] -= pred;
	}

	line[1] -= ( line[2] - line[0] + line[3] ) >> 1;

	for(x=0;x<len;x += 2) {
		A = line[x];
		D = line[x+1];

		a = data[x] = A + ((D+1)>>1);
		data[x+1] = a - D;
	}

}


static void un_tdec_zeros(int *data,int *line,int len)
{
int x,a,A,D;
int E,W,ND,pred;
//int *high;

	/** data is in LLLHHHH order ; we must deshuffle it into line ***/
	/** the HHHH are all zeros, the LLL are goods ***/

	for(x=0;x<(len>>1);x++) {
		line[x+x] = data[x];
	}	

	for(x=len - 2;x>=0;x -= 2) { //** run backwards cuz of ND
		if ( x > 0 ) E = line[x] - line[x - 2];
		else E = 0;
		if ( (x+2)<len ) { W = line[x + 2] - line[x];
			if ( (x+3)<len ) ND = line[x + 3];
			else ND = 0;
		} else W = ND = 0;

		pred = ( E + E + W + W + W + ND + ND ) >> 3;

		line[x+1] = -pred; //** sign difference with encoder
	}

	for(x=0;x<len;x += 2) {
		A = line[x];
		D = line[x+1];

		a = data[x] = A + ((D+1)>>1);
		data[x+1] = a - D;
	}

}

void sp_Transform2D(int **rows, int width, int height, int levels,bool inverse)
{
int x, y, w, h, l;
int *buffer,*tempbuf;

    /* Check the dimensions for compatability. */

    if (width%(1 << (levels+1)) || height%(1 << (levels+1))) {
		errputs("width and height must be divisible by 2^(levels+1)");
		exit(10);
	}
  
    /* Allocate a work array (for transposing columns) */
    
   	if ( (buffer = newarray(int,height+max(width,height))) == NULL ) {
		errputs("malloc failed"); exit(10);
	}
	tempbuf = buffer+height;

    /* do it */
  
	if ( !inverse ) { /* forward transform. */

		for (l = 0; l < levels; l++) {
			w = width >> l;
			h = height >> l;
      
			/* Rows. */
	
			for (y = 0; y < h; y++) {
				do_tdec_line(rows[y],tempbuf,w);
			}
    
			/* Columns. */
	
			for (x = 0; x < w; x++) {
					for (y = 0; y < h; y++) buffer[y] = rows[y][x];
				do_tdec_line(buffer,tempbuf,h);
					for (y = 0; y < h; y++) rows[y][x] = buffer[y];
			}
		}

    } else {

		for (l = levels-1; l >= 0; l--) { /** backwards in scale **/
			w = width >> l;
			h = height >> l;

			/* Columns. */
	
			for (x = 0; x < w; x++) {
					for (y = 0; y < h; y++) buffer[y] = rows[y][x];
				un_tdec_line(buffer,tempbuf,h);
					for (y = 0; y < h; y++) rows[y][x] = buffer[y];
			}

			/* Rows. */
	
			for (y = 0; y < h; y++) {
				un_tdec_line(rows[y],tempbuf,w);
			}    
		}
	}

	free(buffer);
}


void spQuad(int *band,int w,int h,int fullw,bool inverse)
{
int x, y;
int *buffer,*tempbuf,*bptr;

   	if ( (buffer = newarray(int,h+max(w,h))) == NULL ) {
		errputs("malloc failed"); exit(10);
	}
	tempbuf = buffer+h;
  
	if ( !inverse ) { /* forward transform. */
	
		bptr = band;
		for (y = 0; y < h; y++) {
			do_tdec_line(bptr,tempbuf,w);
			bptr += fullw;
		}
    
		for (x = 0; x < w; x++) {
			bptr = band + x;
			for (y = 0; y < h; y++) { buffer[y] = *bptr; bptr += fullw; }
			do_tdec_line(buffer,tempbuf,h);
			bptr = band + x;
			for (y = 0; y < h; y++) { *bptr = buffer[y]; bptr += fullw; }
		}

    } else {

		for (x = 0; x < w; x++) {
			bptr = band + x;
			for (y = 0; y < h; y++) { buffer[y] = *bptr; bptr += fullw; }
			un_tdec_line(buffer,tempbuf,h);
			bptr = band + x;
			for (y = 0; y < h; y++) { *bptr = buffer[y]; bptr += fullw; }
		}

		bptr = band;
		for (y = 0; y < h; y++) {
			un_tdec_line(bptr,tempbuf,w);
			bptr += fullw;
		}

	}

	free(buffer);
}