www.pudn.com > EZW_.rar > WAVELET.C
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Port to C from C++
*
* Mow-Song, Ng 2/9/2002
* msng@mmu.edu.my
* http://www.pesona.mmu.edu.my/~msng
*
* I do not claim copyright to the code, but if you use them or modify them,
* please drop me a mail.
*
*/
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Original copyright info */
/*---------------------------------------------------------------------------*/
// Baseline Wavelet Transform Coder Construction Kit
//
// Geoff Davis
// gdavis@cs.dartmouth.edu
// http://www.cs.dartmouth.edu/~gdavis
//
// Copyright 1996 Geoff Davis 9/11/96
//
// Permission is granted to use this software for research purposes as
// long as this notice stays attached to this software.
//
/*---------------------------------------------------------------------------*/
#include "wavelet.h"
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
WAVELET *WaveletAlloc(FILTERSET *filterset)
{
WAVELET *wavelet;
if((wavelet=(WAVELET *)malloc(sizeof(WAVELET)))==NULL){
return NULL;
}
wavelet->analysisLow=filterset->analysisLow;
wavelet->analysisHigh=filterset->analysisHigh;
wavelet->synthesisLow=filterset->synthesisLow;
wavelet->synthesisHigh=filterset->synthesisHigh;
wavelet->symmetric=filterset->symmetric;
// amount of space to leave for padding vectors for symmetric extensions
wavelet->npad = max(wavelet->analysisLow->size, wavelet->analysisHigh->size);
return wavelet;
}
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletDealloc(WAVELET *wavelet)
{
free(wavelet);
}
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void WaveletTransform1D(WAVELET *wavelet, Real *input, Real *output, int size,
int nsteps, int symExt)
{
int i, currentIndex=0;
Real *data[2];
int lowSize=size, highSize;
// If form of extension unspecified, default to symmetric
// extensions for symmetrical filters and periodic extensions for
// asymmetrical filters
//- Explicitly set symExt = -1
if (symExt == -1){
symExt = wavelet->symmetric;
}
// data[0] and data[1] are padded with npad entries on each end
data [0] = (Real *)calloc(2*wavelet->npad+size, sizeof(Real));
data [1] = (Real *)calloc(2*wavelet->npad+size, sizeof(Real));
for (i = 0; i < size; i++){
data[currentIndex][wavelet->npad+i] = input[i];
}
while(nsteps--){
if (lowSize<=2 && wavelet->symmetric == 1){
WaveletWarning("Reduce # of transform steps or increase signal size.\n");
WaveletWarning("or switch to periodic extension.\n");
WaveletError("Low pass subband is too small\n");
}
// Transform
WaveletTransformStep(wavelet, data[currentIndex], data[1-currentIndex],
lowSize, symExt);
highSize = lowSize/2;
lowSize = (lowSize+1)/2;
// Copy high-pass data to output signal
copy (data[1-currentIndex] + wavelet->npad + lowSize,
output+lowSize, highSize);
// Now pass low-pass data (first 1/2 of signal) back to
// transform routine
currentIndex = 1 - currentIndex;
}
// Copy low-pass data to output signal
copy (data[currentIndex] + wavelet->npad, output, lowSize);
free(data[1]);
free(data[0]);
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletInvert1D(WAVELET *wavelet, Real *input, Real *output, int size,
int nsteps, int symExt)
{
int i;
int currentIndex = 0;
Real *data[2];
int *lowSize, *highSize;
// If form of extension unspecified, default to symmetric
// extensions for symmetrical filters and periodic extensions for
// asymmetrical filters
if (symExt == -1){
symExt = wavelet->symmetric;
}
lowSize = (int *)calloc(nsteps, sizeof(int));
highSize = (int *)calloc(nsteps, sizeof(int));
lowSize[0] = (size+1)/2;
highSize[0] = size/2;
for (i = 1; i < nsteps; i++) {
lowSize[i] = (lowSize[i-1]+1)/2;
highSize[i] = lowSize[i-1]/2;
}
data [0] = (Real*)calloc(2*wavelet->npad+size, sizeof(Real));
data [1] = (Real*)calloc(2*wavelet->npad+size, sizeof(Real));
copy (input, data[currentIndex]+wavelet->npad, lowSize[nsteps-1]);
while (nsteps--) {
// grab the next high-pass component
copy (input + lowSize[nsteps], data[currentIndex]+wavelet->npad+lowSize[nsteps],
highSize[nsteps]);
// Combine low-pass data (first 1/2^n of signal) with high-pass
// data (next 1/2^n of signal) to get higher resolution low-pass data
WaveletInvertStep(wavelet, data[currentIndex], data[1-currentIndex],
lowSize[nsteps]+highSize[nsteps], symExt);
// Now pass low-pass data (first 1/2 of signal) back to
// transform routine
currentIndex = 1 - currentIndex;
}
// Copy inverted signal to output signal
copy (data[currentIndex]+wavelet->npad, output, size);
free(highSize);
free(lowSize);
free(data[1]);
free(data[0]);
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletTransform2D(WAVELET *wavelet, Real *input, Real *output, int hsize, int vsize,
int nsteps, int symExt)
{
int j;
int hLowSize = hsize, hHighSize;
int vLowSize = vsize, vHighSize;
Real *temp_in, *temp_out;
// If form of extension unspecified, default to symmetric
// extensions for symmetrical filters and periodic extensions for
// asymmetrical filters
if (symExt == -1)
symExt = wavelet->symmetric;
temp_in = (Real*)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
temp_out = (Real*)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
copy (input, output, hsize*vsize);
while (nsteps--){
if ((hLowSize <= 2 || vLowSize <= 2) && symExt == 1) {
WaveletWarning ("Reduce # of transform steps or increase signal size.\n");
WaveletWarning ("or switch to periodic extension\n");
WaveletError ("Low pass subband is too small\n");
}
// Do a convolution on the low pass portion of each row
for (j = 0; j < vLowSize; j++) {
// Copy row j to data array
copy (output+(j*hsize), temp_in+wavelet->npad, hLowSize);
// Convolve with low and high pass filters
WaveletTransformStep (wavelet, temp_in, temp_out, hLowSize, symExt);
// Copy back to image
copy (temp_out+wavelet->npad, output+(j*hsize), hLowSize);
}
// Now do a convolution on the low pass portion of each column
for (j = 0; j < hLowSize; j++) {
// Copy column j to data array
copy_p1_skip (output+j, hsize, temp_in+wavelet->npad, vLowSize);
// Convolve with low and high pass filters
WaveletTransformStep (wavelet, temp_in, temp_out, vLowSize, symExt);
// Copy back to image
copy_p2_skip (temp_out+wavelet->npad, output+j, hsize, vLowSize);
}
// Now convolve low-pass portion again
hHighSize = hLowSize/2;
hLowSize = (hLowSize+1)/2;
vHighSize = vLowSize/2;
vLowSize = (vLowSize+1)/2;
}
free(temp_out);
free(temp_in);
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletInvert2D(WAVELET *wavelet, Real *input, Real *output, int hsize, int vsize,
int nsteps, int symExt)
{
int i, j;
int *hLowSize, *hHighSize, *vLowSize, *vHighSize;
Real *temp_in, *temp_out;
// If form of extension unspecified, default to symmetric
// extensions for symmetrical filters and periodic extensions for
// asymmetrical filters
if (symExt == -1)
symExt = wavelet->symmetric;
hLowSize = (int *) calloc(nsteps, sizeof(int));
hHighSize = (int *) calloc(nsteps, sizeof(int));
vLowSize =(int *) calloc(nsteps, sizeof(int));
vHighSize =(int *) calloc(nsteps, sizeof(int));
hLowSize[0] = (hsize+1)/2;
hHighSize[0] = hsize/2;
vLowSize[0] = (vsize+1)/2;
vHighSize[0] = vsize/2;
for (i = 1; i < nsteps; i++) {
hLowSize[i] = (hLowSize[i-1]+1)/2;
hHighSize[i] = hLowSize[i-1]/2;
vLowSize[i] = (vLowSize[i-1]+1)/2;
vHighSize[i] = vLowSize[i-1]/2;
}
temp_in = (Real *)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
temp_out = (Real *)calloc(2*wavelet->npad+MAX(hsize,vsize), sizeof(Real));
copy (input, output, hsize*vsize);
while (nsteps--) {
// Do a reconstruction for each of the columns
for (j = 0; j < hLowSize[nsteps]+hHighSize[nsteps]; j++) {
// Copy column j to data array
copy_p1_skip (output+j, hsize, temp_in+wavelet->npad,
vLowSize[nsteps]+vHighSize[nsteps]);
// Combine low-pass data (first 1/2^n of signal) with high-pass
// data (next 1/2^n of signal) to get higher resolution low-pass data
WaveletInvertStep (wavelet, temp_in, temp_out,
vLowSize[nsteps]+vHighSize[nsteps], symExt);
// Copy back to image
copy_p2_skip (temp_out+wavelet->npad, output+j, hsize,
vLowSize[nsteps]+vHighSize[nsteps]);
}
// Now do a reconstruction pass for each row
for (j = 0; j < vLowSize[nsteps]+vHighSize[nsteps]; j++) {
// Copy row j to data array
copy (output + (j*hsize), temp_in+wavelet->npad,
hLowSize[nsteps]+hHighSize[nsteps]);
// Combine low-pass data (first 1/2^n of signal) with high-pass
// data (next 1/2^n of signal) to get higher resolution low-pass data
WaveletInvertStep (wavelet, temp_in, temp_out,
hLowSize[nsteps]+hHighSize[nsteps], symExt);
// Copy back to image
copy (temp_out+wavelet->npad, output + (j*hsize),
hLowSize[nsteps]+hHighSize[nsteps]);
}
}
free(hLowSize);
free(hHighSize);
free(vLowSize);
free(vHighSize);
free(temp_in);
free(temp_out);
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
// Do symmetric extension of data using prescribed symmetries
// Original values are in output[npad] through output[npad+size-1]
// New values will be placed in output[0] through output[npad] and in
// output[npad+size] through output[2*npad+size-1] (note: end values may
// not be filled in)
// left_ext = 1 -> extension at left bdry is ... 3 2 1 | 0 1 2 3 ...
// left_ext = 2 -> extension at left bdry is ... 3 2 1 0 | 0 1 2 3 ...
// right_ext = 1 or 2 has similar effects at the right boundary
//
// symmetry = 1 -> extend symmetrically
// symmetry = -1 -> extend antisymmetrically
void WaveletSymmetricExtension(WAVELET *wavelet, Real *output, int size, int leftExt, int
rightExt, int symmetry)
{
int i;
int first = wavelet->npad, last = wavelet->npad + size-1;
int originalFirst, originalLast, originalSize, period;
int nextend;
if (symmetry == -1) {
if (leftExt == 1)
output[--first] = 0;
if (rightExt == 1)
output[++last] = 0;
}
originalFirst = first;
originalLast = last;
originalSize = originalLast-originalFirst+1;
period = 2 * (last - first + 1) - (leftExt == 1) - (rightExt == 1);
if (leftExt == 2)
output[--first] = symmetry*output[originalFirst];
if (rightExt == 2)
output[++last] = symmetry*output[originalLast];
// extend left end
nextend = MIN (originalSize-2, first);
for (i = 0; i < nextend; i++) {
output[--first] = symmetry*output[originalFirst+1+i];
}
// should have full period now -- extend periodically
while (first > 0) {
first--;
output[first] = output[first+period];
}
// extend right end
nextend = MIN (originalSize-2, 2*wavelet->npad+size-1 - last);
for (i = 0; i < nextend; i++) {
output[++last] = symmetry*output[originalLast-1-i];
}
// should have full period now -- extend periodically
while (last < 2*wavelet->npad+size-1) {
last++;
output[last] = output[last-period];
}
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
// Do periodic extension of data using prescribed symmetries
// Original values are in output[npad] through output[npad+size-1]
// New values will be placed in output[0] through output[npad] and in
// output[npad+size] through output[2*npad+size-1] (note: end values may
// not be filled in)
void WaveletPeriodicExtension(WAVELET *wavelet, Real *output, int size)
{
int first = wavelet->npad, last = wavelet->npad + size-1;
// extend left periodically
while (first > 0) {
first--;
output[first] = output[first+size];
}
// extend right periodically
while (last < 2*wavelet->npad+size-1) {
last++;
output[last] = output[last-size];
}
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletTransformStep (WAVELET *wavelet, Real *input, Real *output,
int size, int symExt)
{
int i, j;
int lowSize=(size+1)/2;
int leftExt, rightExt;
if (wavelet->analysisLow->size%2){
// odd filter length
leftExt=rightExt=1;
}
else{
leftExt=rightExt=2;
}
if (symExt){
WaveletSymmetricExtension(wavelet, input, size, leftExt, rightExt, 1);
}
else{
WaveletPeriodicExtension(wavelet, input, size);
}
// coarse detail
// xxxxxxxxxxxxxxxx --> HHHHHHHHGGGGGGGG
for (i = 0; i < lowSize; i++) {
output[wavelet->npad+i] = 0.0;
for (j = 0; j < wavelet->analysisLow->size; j++) {
output [wavelet->npad+i] +=
input[wavelet->npad + 2*i + wavelet->analysisLow->firstIndex + j] *
wavelet->analysisLow->coeff[j];
}
}
for (i = lowSize; i < size; i++) {
output[wavelet->npad+i] = 0.0;
for (j = 0; j < wavelet->analysisHigh->size; j++) {
output [wavelet->npad+i] +=
input[wavelet->npad + 2*(i-lowSize) + wavelet->analysisHigh->firstIndex + j] *
wavelet->analysisHigh->coeff[j];
}
}
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletInvertStep (WAVELET *wavelet, Real *input, Real *output,
int size, int symExt)
{
int i, j;
int leftExt, rightExt, symmetry;
Real *temp;
int firstIndex, lastIndex;
// amount of low and high pass -- if odd # of values, extra will be
// low pass
int lowSize = (size+1)/2, highSize = size/2;
symmetry = 1;
if (wavelet->analysisLow->size % 2 == 0) {
// even length filter -- do (2, X) extension
leftExt = 2;
} else {
// odd length filter -- do (1, X) extension
leftExt = 1;
}
if (size % 2 == 0) {
// even length signal -- do (X, 2) extension
rightExt = 2;
} else {
// odd length signal -- do (X, 1) extension
rightExt = 1;
}
temp = (Real *)calloc(2*wavelet->npad+lowSize, sizeof(Real));
for (i = 0; i < lowSize; i++) {
temp[wavelet->npad+i] = input[wavelet->npad+i];
}
if (symExt){
WaveletSymmetricExtension (wavelet, temp, lowSize, leftExt, rightExt, symmetry);
}
else{
WaveletPeriodicExtension (wavelet, temp, lowSize);
}
// coarse detail
// HHHHHHHHGGGGGGGG --> xxxxxxxxxxxxxxxx
for (i = 0; i < 2*wavelet->npad+size; i++){
output[i] = 0.0;
}
firstIndex = wavelet->synthesisLow->firstIndex;
lastIndex = wavelet->synthesisLow->size - 1 + firstIndex;
for (i = -lastIndex/2; i <= (size-1-firstIndex)/2; i++) {
for (j = 0; j < wavelet->synthesisLow->size; j++) {
output[wavelet->npad + 2*i + firstIndex + j] +=
temp[wavelet->npad+i] * wavelet->synthesisLow->coeff[j];
}
}
leftExt = 2;
if (wavelet->analysisLow->size % 2 == 0) {
// even length filters
rightExt = (size % 2 == 0) ? 2 : 1;
symmetry = -1;
}
else {
// odd length filters
rightExt = (size % 2 == 0) ? 1 : 2;
symmetry = 1;
}
for (i = 0; i < highSize; i++) {
temp[wavelet->npad+i] = input[wavelet->npad+lowSize+i];
}
if (symExt){
WaveletSymmetricExtension (wavelet, temp, highSize, leftExt, rightExt,
symmetry);
}
else{
WaveletPeriodicExtension (wavelet, temp, highSize);
}
firstIndex = wavelet->synthesisHigh->firstIndex;
lastIndex = wavelet->synthesisHigh->size - 1 + firstIndex;
for (i = -lastIndex/2; i <= (size-1-firstIndex)/2; i++) {
for (j = 0; j < wavelet->synthesisHigh->size; j++) {
output[wavelet->npad + 2*i + firstIndex + j] +=
temp[wavelet->npad+i] * wavelet->synthesisHigh->coeff[j];
}
}
free(temp);
}
// copy length elements from p1 to p2
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void copy (const Real *p1, Real *p2, const int length){
int temp = length;
while(temp--) *p2++ = *p1++;
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void copy_p1_skip (const Real *p1, const int stride1, Real *p2, const int length){
int temp = length;
while(temp--){
*p2++ = *p1;
p1 += stride1;
}
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void copy_p2_skip (const Real *p1, Real *p2, const int stride2, const int length){
int temp = length;
while(temp--){
*p2 = *p1++;
p2 += stride2;
}
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletWarning(char *fmt, ...)
{
va_list argptr;
va_start( argptr, fmt );
fprintf( stderr, "WaveletWarning: " );
vprintf( fmt, argptr );
va_end( argptr );
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void WaveletError(char *fmt, ...)
{
va_list argptr;
va_start( argptr, fmt );
fprintf(stderr, "WaveletError: " );
vprintf( fmt, argptr );
va_end( argptr );
exit( -1 );
}