wavelet.cpp

wavlet compression on c++ only cods

// Copyright (c) 1999 Stephen T. Welstead. All rights reserved.
// wavelet.cpp Wavelet classes.

#include <stdio.h>
#include <math.h>

#include "utofile.h"

#include "utm.h"
#include "wavelet.h"
#include "waveproc.h"

// Assumes all data arrays start with index 1 in each dimension.
// Length of all data arrays must be a power of 2.

void twavelet::transform (float *data_vector,int n) {
// Note: n must be a power of 2!
     int nn;
     if (n < 4) return;
     for (nn = n; nn >= 4; nn >>= 1)
	 filter (data_vector,nn);
     return;
     }

void twavelet::inverse_transform (float *data_vector,int n) {
// Note: n must be a power of 2!
     int nn;
     if (n < 4) return;
     for (nn = 4; nn <= n; nn <<= 1)
	 transpose_filter (data_vector,nn);
     return;
     }

int twavelet::filter(float *data,int no_of_pts) {
     // Assumes data allocated 1..no_of_pts and
     // coeff allocated 0..no_of_coeffs-1
     int i,j,index_l,index_h;
     int half_no_of_pts = no_of_pts/2;
     int sign;
     double sum_l,sum_h;
     double *x;
     if (no_of_pts < no_of_coeffs) return 0;
     if (!(x=allocate_d_vector(1,no_of_pts))) return 0;
     for (i=1;i<=half_no_of_pts;i++) {
	 sum_l = 0.0;
	 sum_h = 0.0;
	 sign = 1;
	 for (j=1;j<=no_of_pts;j++) {
	     index_l = j-2*i+1;
	     index_h = 2*i-j;
	     while (index_l < 0)
		   index_l += no_of_pts;
	     if (index_l < no_of_coeffs)
		sum_l += coeff[index_l]*data[j];
	     while (index_h < 0)
		   index_h += no_of_pts;
	     if (index_h < no_of_coeffs)
		sum_h += sign*coeff[index_h]*data[j];
	     sign = -sign;
	     }  // end j
	x[i] = sum_l;
	x[i+half_no_of_pts] = sum_h;
	} // end i
     for (i=1;i<=no_of_pts;i++)
	 data[i] = x[i];
     free_d_vector(x,1);
     return 1;
     }

int twavelet::transpose_filter (float *data,int no_of_pts) {
     // Assumes data allocated 1..no_of_pts and
     // coeff allocated 0..no_of_coeffs-1
     int i,j,index_l,index_h;
     int half_no_of_pts = no_of_pts/2;
     int sign;
     double sum_l,sum_h;
     double *x;
     if (no_of_pts < no_of_coeffs) return 0;
     if (!(x=allocate_d_vector(1,no_of_pts))) return 0;
     sign = 1;
     for (i=1;i<=no_of_pts;i++) {
	 sum_l = 0.0;
	 sum_h = 0.0;
	 for (j=1;j<=half_no_of_pts;j++) {
	     index_l = i-2*j+1;
	     index_h = 2*j-i /* -2+no_of_coeffs */ ;
	     while (index_l < 0)
		   index_l += no_of_pts;
	     while (index_h < 0)
		   index_h += no_of_pts;
	     if (index_l < no_of_coeffs)
		sum_l += coeff[index_l]*data[j];
	     if (index_h < no_of_coeffs)
		sum_h += sign*coeff[index_h]*data[j+half_no_of_pts];
	     }  // end j
	sign = -sign;
	x[i] = sum_l + sum_h;
	} // end i
     for (i=1;i<=no_of_pts;i++)
	 data[i] = x[i];
     free_d_vector(x,1);
     return 1;
     }

double sqrt_2_inv = 1.0/sqrt(2.0);
double Haar_coeff[] = {sqrt_2_inv,sqrt_2_inv};

tHaar_wavelet::tHaar_wavelet(void):twavelet (Haar_coeff,2){};

void tHaar_wavelet::transform (float *data_vector,int n) {
// Note: n must be a power of 2!
     int i;
     i = n;
     while (i >=2) {
	 filter (data_vector,i);
	 i /= 2;
	 }
     return;
     }

void tHaar_wavelet::inverse_transform (float *data_vector,int n) {
// Note: n must be a power of 2!
     int i;
     i = 2;
     while (i <= n) {
	 transpose_filter (data_vector,i);
	 i *= 2;
	 }
     return;
     }

double D4_coeff[] = {D4_0,D4_1,D4_2,D4_3};
tDaub4_wavelet::tDaub4_wavelet(void):twavelet (D4_coeff,4){};

double D6_coeff[] = {D6_0,D6_1,D6_2,D6_3,D6_4,D6_5};
tDaub6_wavelet::tDaub6_wavelet(void):twavelet (D6_coeff,6){};

twavelet_2d_array::twavelet_2d_array (float **the_values,
	 twavelet *the_wavelet,int the_rows,int the_cols) {
     wavelet = the_wavelet;
     values = the_values;
     rows = the_rows;
     cols = the_cols;
     }

int write_coeffs (char *filename,float **the_values,
      int print_rows,int print_cols) {
    FILE *outfile;
    int i,j;
    if (!open_output_text_file (&outfile,filename)) return 0;
    for (i=1;i<=print_rows;i++) {
	for (j=1;j<=print_cols;j++)
	    fprintf (outfile,"%7.2f  ",the_values[i][j]);
	fprintf (outfile,"\n");
	}
    fclose (outfile);
    return 1;
    }

void twavelet_2d_array::transform () {
    float *col_values;
    for (int i= /*2*/ 1;i<=rows;i++) {
	wavelet->transform (values[i],cols);
	}  // end i
    col_values = allocate_f_vector (1,rows);
    for (int j=1;j<=cols;j++) {
	for (i=1;i<=rows;i++)
	    col_values[i] = values[i][j];
	wavelet->transform (col_values,rows);
	for (i=1;i<=rows;i++)
	    values[i][j] = col_values[i];
	} // end j
    free_f_vector (col_values,1);
    lowpass_value = values[1][1];
    transform_max = transform_min = values[1][2];
    for (i= 1;i<=rows;i++)
	for (j=1;j<=cols;j++) {
            // Don't include lowpass_value in max-min comparisons:
	    if ((i!=1)||(j!=1)) {
	    if (values[i][j] > transform_max) transform_max = values[i][j];
	    if (values[i][j] < transform_min) transform_min = values[i][j];
	    }
	    }  // end i,j
    return;
    }

void twavelet_2d_array::inverse_transform () {
    float *col_values;
    int i,j;
    col_values = allocate_f_vector (1,rows);
    for (j=1;j<=cols;j++) {
	for (i=1;i<=rows;i++)
	    col_values[i] = values[i][j];
	wavelet->inverse_transform (col_values,rows);
	for (i=1;i<=rows;i++)
	    values[i][j] = col_values[i];
	} // end j
    free_f_vector (col_values,1);
    for (i= /*2*/ 1;i<=rows;i++) {
	wavelet->inverse_transform (values[i],cols);
	}  // end i
    inverse_max = inverse_min = values[1][1];
    for (i=1;i<=rows;i++)
	for (j=1;j<=cols;j++) {
	    if (values[i][j] > inverse_max) inverse_max = values[i][j];
	    if (values[i][j] < inverse_min) inverse_min = values[i][j];
	    }  // end i,j
    return;
    }

void twavelet_2d_array::decimate (float percent) {
    // Set all but top 'percent' array values equal to 0.
    unsigned long n,nrows = rows,ncols = cols;
    if (percent >= 99.99) return;  
    n = (percent/100.0)*nrows*ncols;
    decimate_array (values,rows,cols,n);
    return;
    }