/* * * GWIC * * (c) Joonas Lehtinen (jole@jole.fi), 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);}