lifting53.cpp

双正交5/3小波的提升格式构造

/*
  this program performs 1D and 2D wavelet decomposition 
  and synthesis using lifting algorithm. the wavelet used
  here is a set of invertible wavelet that map integer to
  integer, namely 5/3 wavelet.

  forward transformation:
  d[n] = d_0[n] - \leftfloor 1/2 (s_0[n+1] + s_0[n]) \rightfloor
  s[n] = s_0[n] + \leftfloor 1/4 (d[n] + d[n-1]) + 1/2 \rightfloor

  written by Y.B. Mao at Dept. of Automation, NJUST
  Oct. 23, 2005
  Version 1.0

  All rights reserved.  
*/
# include <stdlib.h>
# include <stdio.h>
# include <malloc.h>
# include "lifting53.h"

/*
 one round 5/3 wavelet decomposition using lifting
*/
int wavelet_analysis( int * data, int n,   // original data
				      int * approximation, // approximation coefficients
				      int * detail )       // detail coefficients
{
	int m, i;
	int * x, * d;
	
	m = n/2;
	n = m * 2;  /* n must be even */

	x = ( int * ) malloc( sizeof(int) * (n+3) );
	d = ( int * ) malloc( sizeof(int) * (m+1) );

	for ( i = 0; i < n; i++ ) /* boundary extension, CB ABCDEF E */
	{
		x[i+2] = data[i];
	}
	x[0] = data[2];
	x[1] = data[1];
	x[n+2] = data[n-2];
	/* detail (odd) */
	for ( i = 0; i <= m; i++ )
	{
		d[i] = x[i*2+1] - (x[i*2]+x[i*2+2])/2;
	}
	for ( i = 0; i < m; i++ )
		detail[i] = d[i+1];
	/* approximation (even) */
	for ( i = 0; i < m; i++ )
	{
		approximation[i] = x[2*i+2] + ( d[i+1]+d[i]+2 )/4;
	}

	free( x );
	free( d );

	return( m );
}

/*
 one round wavelet synthesis
*/
int wavelet_synthesis( int * approximation, // approximation coefficients
				       int * detail,        // detail coefficients
					   int n,               // n is length of both detail and approximation
				       int * data )         // synthesised data
{
	int * s, * d, * x;
	int i, m;

	m = n * 2;
	s = (int *) malloc( sizeof(int) * (n+1) );
	d = (int *) malloc( sizeof(int) * (n+2) );
	x = (int *) malloc( sizeof(int) * (m+1) );

	/* boundary extension */
	for ( i = 0; i < n; i++ ) s[i] = approximation[i];
	s[n] = approximation[n-1]; /* right extension, ...EF F */
	for ( i = 0; i < n; i++ ) d[i+1] = detail[i];
	d[0] = detail[0];  /* left extension, A ABC ... */
	d[n+1] = detail[n-2]; /* right extension, ... DEF E */
	/* lifting */
	for ( i = 0; i < n+1; i++ )
		x[i*2] = s[i]-(d[i]+d[i+1]+2)/4;
	for ( i = 0; i < n; i++ )
		x[i*2+1] = d[i+1] + (x[i*2]+x[i*2+2])/2;
	/* copy x to data */
	for ( i = 0; i < m; i++ )
		data[i] = x[i];

	free( x );
	free( s );
	free( d );

	return( m );
}

/*
 multi-level decomposition for 1D signal
*/
int multilevel_decomposition_1D( int * original, // original signal
								 int n,             // number of pointers, should be multiple of 2
								 int * resultdata,    // decomposited data
								 int levels )       // decomposition levels  
{
	int i, x, nn;
	int rtv;
	int * l, * h, * data;

	nn = (n>>levels)<<levels; // guartee nn to be multiple of 2
	for ( x = nn; x < n; x++ ) resultdata[x] = 0;
	l = new int[n];
	h = new int[n];
	data = new int[n];
	for ( x = 0; x < nn; x++ ) l[x] = original[x];
	for ( i = 0; i < levels; i++ )
	{
		for ( x = 0; x < nn; x++ ) data[x] = l[x];
		rtv = wavelet_analysis( data, nn, l, h );
		nn = nn/2;
		for ( x = 0; x < nn; x++ ) resultdata[x] = l[x];
		for ( x = 0; x < nn; x++ ) resultdata[nn+x] = h[x];
	}
	delete[] l;
	delete[] h;
	delete[] data;

	return( nn ); // return total points subject to decomposition
}

/*
 multi-level reconstruction for 1D signal
*/
int multilevel_reconstruction_1D( int * de_data, // deconstructed signal
								  int n,             // number of pointers, should be multiple of 2
								  int * resultdata,    // reconstructed data
								  int levels )       // reconstruct levels  
{
	int i, x, nn;
	int rtv;
	int * l, * h, * data;

	nn = n>>levels;
	l = new int[n];
	h = new int[n];
	data = new int[n];
	for ( x = 0; x < n; x++ ) data[x] = de_data[x];
	for ( i = 0; i < levels; i++ )
	{
		for ( x = 0; x < nn; x++ ) l[x] = data[x];
		for ( x = 0; x < nn; x++ ) h[x] = data[nn+x];
		rtv = wavelet_synthesis( l, h, nn, data );
		nn = nn*2;
	}
	for ( x = 0; x < n; x++ ) resultdata[x] = data[x];
	delete[] l;
	delete[] h;
	delete[] data;

	return( nn ); // return total points subject to decomposition
}

/*
 multi-level decomposition for 2D signal
*/
int multilevel_decomposition_2D( int * image,         // original signal
								 int width, int height,  // number of pointers, should be multiple of 2
								 int * wavelet_data,  // decomposited data
								 int levels )            // decomposition levels  
{
	int i, x, y, ww, hh;
	int rtv;
	int * l, * h, * data;

	ww = ( width>>levels )<<levels; // guartee nn to be multiple of 2
	hh = ( height>>levels )<<levels;

	for ( y = 0; y < height; y++ )
		for ( x = 0; x < width; x++ )
			wavelet_data[y*width+x] = image[y*width+x];

	for ( i = 0; i < levels; i++ )
	{
		for ( y = 0; y < hh; y++ )
		{
			l = new int[ww];
			h = new int[ww];
			data = new int[ww];
			for ( x = 0; x < ww; x++ ) data[x] = wavelet_data[ y*width+x ];
			rtv = wavelet_analysis( data, ww, l, h );
			for ( x = 0; x < ww/2; x++ ) wavelet_data[y*width+x] = l[x];
			for ( x = 0; x < ww/2; x++ ) wavelet_data[y*width+x+ww/2] = h[x];
			delete[] l;
			delete[] h;
			delete[] data;
		}
		for ( x = 0; x < ww; x++ )
		{
			l = new int[hh];
			h = new int[hh];
			data = new int[hh];
			for ( y = 0; y < hh; y++ ) data[y] = wavelet_data[ y*width+x ];
			rtv = wavelet_analysis( data, hh, l, h );
			for ( y = 0; y < hh/2; y++ ) wavelet_data[y*width+x] = l[y];
			for ( y = 0; y < hh/2; y++ ) wavelet_data[(y+hh/2)*width+x] = h[y];
			delete[] l;
			delete[] h;
			delete[] data;
		}
		ww = ww/2;
		hh = hh/2;
	}

	return( levels ); // return decomposition levels
}

/*
 multi-level reconstruction for 2D signal
*/
int multilevel_reconstruction_2D( int * wavelet_data,  // wavelet coefficients
								  int width, int height,  // number of pointers, should be multiple of 2
								  int * image,         // reconstructed data
								  int levels )            // reconstruction levels
{
	int i, x, y, ww, hh;
	int rtv;
	int * l, * h, * data;

	ww = width>>levels;  // guartee nn to be multiple of 2
	hh = height>>levels;

	for ( y = 0; y < height; y++ )
		for ( x = 0; x < width; x++ )
			image[y*width+x] = wavelet_data[y*width+x];

	for ( i = 0; i < levels; i++ )
	{

		for ( x = 0; x < ww*2; x++ )
		{
			l = new int[hh];
			h = new int[hh];
			data = new int[hh*2];
			for ( y = 0; y < hh; y++ ) 
			{
				l[y] = image[ y*width+x ];
				h[y] = image[ (y+hh)*width+x ];
			}
			rtv = wavelet_synthesis( l, h, hh, data );
			for ( y = 0; y < hh*2; y++ ) image[y*width+x] = data[y];
			delete[] l;
			delete[] h;
			delete[] data;
		}
		for ( y = 0; y < hh*2; y++ )
		{
			l = new int[ww];
			h = new int[ww];
			data = new int[ww*2];
			for ( x = 0; x < ww; x++ )
			{
				l[x] = image[ y*width+x ];
				h[x] = image[ y*width+x+ww ];
			}
			rtv = wavelet_synthesis( l, h, ww, data );
			for ( x = 0; x < ww*2; x++ ) image[y*width+x] = data[x];
			delete[] l;
			delete[] h;
			delete[] data;
		}
		ww = ww*2;
		hh = hh*2;
	}

	return( levels ); // return decomposition levels
}