transform.c

一种新的机遇小波变换的压缩存储方法

/*
 *
 * GWIC
 *
 * (c) Joonas Lehtinen ([e-mail]), TUCS, 1998
 *
 */

#include "gwic.h"

/* Daubechies filters, N = 2 (Vetterli: Subband.., p.260) */

const float tr_daub2_high[4] = {
  .129409522, .22414386, -.8365163, .48296291 };
const float tr_daub2_low[4] = { 
  .48296291 , .8365163, .22414386, -.129409522 };

/* Biorthogonal B97 filters */
const float tr_b97_f_high[6] = { 
  0.788485, -0.418092, -0.040690,  0.064539, 0.0, 0.0 };
const float tr_b97_i_high[6] = { 
  0.852699,  -0.377403, -0.110624, 0.023849, 0.037829, 0.0 };
const float tr_b97_f_low[6] = { 
  0.852699,  0.377403, -0.110624, -0.023849, 0.037829, 0.0 };
const float tr_b97_i_low[6] = { 
  0.788485, 0.418092, -0.040690,  -0.064539, 0.0, 0.0 };

/* Filter float vector in as defined by transform */
void filter_vector(float *in, float *out, int l, BYTE transform)
{
  int i,j,k,h,n,m;
  
  switch(transform) {
  case TR_DAUB2:
    h = l >> 1;
    for (i=0; i<h; i++) {
      out[i] = out[i+h] = 0;
      for (j=0; j<4; j++) {
	k = (i << 1) + j; if (k>=l) k-=l; 
	out[i] += in[k] * tr_daub2_low[j];
	out[i+h] += in[k] * tr_daub2_high[j];
      }
    }
    break;
  case TR_B97:
    h = l >> 1;
    for (i=0; i<h; i++) {
      k = i<<1;
      out[i] = in[k] * tr_b97_f_low[0];
      out[i+h] = in[k+1] * tr_b97_f_high[0];
      for (j=1; j<5; j++) {
	n = k + j; if (n>=l) n-=l; 
	m = k - j; if (m<0) m+=l; 
	out[i] += (in[n]+in[m]) * tr_b97_f_low[j];
	n = k + j +1; if (n>=l) n-=l; 
	m = k - j +1; if (m<0) m+=l; 
	out[i+h] += (in[n]+in[m])* tr_b97_f_high[j];
      }
    }
    break;
  default: 
    TERMINATE("Unsupported transform requested");
  }
}

/* Inverse-filter float vector in as defined by transform */
void inverse_filter_vector(float *in, float *out, int l, BYTE transform)
{
  int i,j,k,h,n,m;

  switch(transform) {
  case TR_DAUB2:
    h=l >> 1;
    for(i=0;i<l;i++) out[i]=0.0;
    for (i=0; i<h; i++) {
      for (j=0; j<4; j++) {
	k = (i<<1)+j; if (k>=l) k-=l; 
	out[k] += in[i]*tr_daub2_low[j] + in[i+h]*tr_daub2_high[j];
      }
    }
    break;
  case TR_B97:
    h=l>>1;
    for(i=0;i<l;i++) out[i]=0.0;
    for (i=0; i<h; i++) {
      k = i<<1;
      out[k] += in[i] * tr_b97_i_low[0];
      out[k+1] += in[i+h] * tr_b97_i_high[0];
      for (j=1; j<5; j++) {
	n = k + j; if (n>=l) n-=l; 
	m = k - j; if (m<0) m+=l; 
	out[n] += in[i] * tr_b97_i_low[j];
	out[m] += in[i] * tr_b97_i_low[j];
	n = k + j +1; if (n>=l) n-=l; 
	m = k - j +1; if (m<0) m+=l; 
	out[n] += in[i+h] * tr_b97_i_high[j];
	out[m] += in[i+h] * tr_b97_i_high[j];
      }
    }
    break;
  default: 
    TERMINATE("Unsupported transform requested");
  }
}

/*D*/
void debug_b97() 
{
  float *b, *t, *e, *ib, *ob;
  int i,j;

  MALLOC(ib,100*sizeof(float));
  MALLOC(ob,100*sizeof(float));
  ib+=20;
  ob+=20;
  MALLOC(b,100*sizeof(float));
  MALLOC(t,100*sizeof(float));
  MALLOC(e,100*sizeof(float));
  
  for (i=0;i<20;i++) *(ib+i) = *(b+i) = 10+i+2*(i&1);
  *(ib+10) = *(b+10) = 20;
  filter_vector(ib,ob,20,TR_B97);
  for (i=0;i<20;i++) *(t+i) = *(ib+i) = *(ob+i);
  inverse_filter_vector(ib,ob,20,TR_B97);
  for (i=0;i<20;i++) *(e+i) = *(ob+i); 
  

  for(i=0; i<20; i++) {
    printf("%i:\t%f\t%f\t%f\n",i,*(b+i),*(t+i),*(e+i));
  }
  
}

int forward_transform(float *table, int width, int height, BYTE transform)
{
  int w,h,i,j;
  float *ibuf,*obuf;
  float *p;
  int levels;

  w = width; h = height; levels=0;
  MALLOC(ibuf,((w>h?w:h)+20)*sizeof(float));
  MALLOC(obuf,((w>h?w:h)+20)*sizeof(float));
  ibuf += 10; obuf += 10;

  /* levels */
  while(!((w|h)&1) && w>=8 && h>=8) { 
    /* rows */
    for(j=0; j<h; j++) {
      for(p=table+j*width, i=0; i<w;i++) *(ibuf+i) = *p++;
      filter_vector(ibuf,obuf,w,transform);
      for(p=table+j*width, i=0; i<w;i++) *p++ = *(obuf+i);
    }
    /* cols */
    for(i=0; i<w; i++) {
      for(p=table+i, j=0; j<h;j++,p+=width) *(ibuf+j) = *p;
      filter_vector(ibuf,obuf,h,transform);
      for(p=table+i, j=0; j<h;j++,p+=width) *p = *(obuf+j);
    }
    w = w >> 1 ; h = h >> 1; levels++;
  }

  ibuf -= 10; obuf -= 10;
  free(ibuf); free(obuf);

  return levels;
}

void inverse_transform(float *table, int width, int height, BYTE transform)
{
  int w,h,i,j;
  float *ibuf,*obuf;
  float *p;

  w = width; h = height;
  MALLOC(ibuf,((w>h?w:h)+20)*sizeof(float));
  MALLOC(obuf,((w>h?w:h)+20)*sizeof(float));
  ibuf += 10; obuf += 10;

  while(!(((w>>1)|(h>>1))&1) && (w>>1)>=8 && (h>>1)>=8) { w=w>>1; h=h>>1; }
  
  /* levels */
  while(w<=width && h<=height) {
    /* cols */
    for(i=0; i<w; i++) {
      for(p=table+i, j=0; j<h;j++,p+=width) *(ibuf+j) = *p;
      inverse_filter_vector(ibuf,obuf,h,transform);
      for(p=table+i, j=0; j<h;j++,p+=width) *p = *(obuf+j);
    }

    /* rows */
    for(j=0; j<h; j++) {
      for(p=table+j*width, i=0; i<w;i++) *(ibuf+i) = *p++;
      inverse_filter_vector(ibuf,obuf,w,transform);
      for(p=table+j*width, i=0; i<w;i++) *p++ = *(obuf+i);
    }

    w = w << 1 ; h = h << 1; 
  }

  ibuf -= 10; obuf -= 10;
  free(ibuf); free(obuf);
}