Link to home
Start Free TrialLog in
Avatar of AdamB_Open
AdamB_Open

asked on

Wavelets

I need an algorith to do a forward wavelet transformation of 2-D data using Gabor wavelets.  The algorithm will be implemented in a specialised language, but I can read VB,C/C++.  I can find plenty of mathematical explanations and pseudocode for wavelets on the net, but I have difficulty converting this into code.
Avatar of cookre
cookre
Flag of United States of America image
If you don't get any direct responses, you may want to post some of the pseudo code for folks here to work on.
Avatar of atlan2
atlan2
"In image analysis, the two-dimensional wavelets are called the Gabor functions, and the correspondig translationally invariant features are the Gabor amplitude transforms" (from T.Kohonen book: Self-Organizing Maps. 2.nd ed. p.178, Springer)


I haven't tested this source code yet, but it's short enough to paste here.


/*************************************************************
 *      Title : wavelet.c
 *
 *  Basic set of subroutines for 1D and 2D wavelet transformation.
 *
 *  Michael Unser / BEIP      OCT-93
 *
 *************************************************************/
#include <math.h>
#include <stdlib.h>
#define MAX(x,y) ( ( (x) > (y) ) ? (x) : (y) )
#define MIN(x,y) ( ( (x) < (y) ) ? (x) : (y) )
#define ROUND( a )      ( ( (a) < 0 ) ? (int) ( (a) - 0.5 ) : \
/*************************************************************
 *  FUNCTION : wavelet_transform
 *  PURPOSE : Performs a wavelet transform with a depth of abs(nlevels)
 *  using the lowpass and highpass filters v[] and w[].
 *  nlevels > 0 : direct transform
 *  nlevels < 0 : indirect transform
 *
 *  Michael Unser / BEIP      OCT-93
 *                                          OCT-94  (mean is added back to lowpass channel)
 *                                          OCT-94  (added periodic boundary conditions -> mirror==0)
 *
 *************************************************************/
  void ALRT(char      *str);
  void copyr(float *ix1,int nx1,int ny1, float *ix2,int nx2,int ny2);
  void uqmfup(float y1[],float y2[],int ny,float x[],int *nx,float v[],int
nv,float w[],int nw);
  void uqmfdown(float x[],int nx,float y1[],float y2[],int *ny,float
v[],int nv,float w[],int nw);
  void uqmfupp(float y1[],float y2[],int ny,float x[],int *nx,float v[],int
nv,float w[],int nw);
  void uqmfdownp(float x[],int nx,float y1[],float y2[],int *ny,float
v[],int nv,float w[],int nw);

int wavelet_transform(float *ix,float *wave,int M,int N,float v[],int
nv,float w[],int nw,int nlevels,float *mean,int mirror);
int wavelet_transform(float *ix,float *wave,int M,int N,float v[],int
nv,float w[],int nw,int nlevels,float *mean,int mirror)
{

  long maxsize,kk,kmax;
  int nx,ny;
  short i,j,M1,N1,M2,N2,k;
  float *p1,*p2,*p3,*p11,*p22,*p33;
  float *pvv,*pvw,*pwv,*pww;
  float *x,*y1=0,*y2=0;
  float xm=0;
  short levels=abs(nlevels),reduc=pow(2.,(double)levels);

  if (nlevels==0) {
        copyr(ix,M,N,wave,M,N);
        return 1;
  }

  M1=(M/reduc)*reduc;
  N1=(N/reduc)*reduc;

  if ((M1!=M) || (N1!=N)) {
        ALRT("Non-compatible dimension of input image for wavelet transform");
        return (0);
  }
  if ((M1<2) || (N1<2)) {
        ALRT("Can't handle so many levels");
        return (0);
  }

  maxsize=MAX(M,N)*sizeof(float);
  x=(float*)malloc((size_t)maxsize);
  y1=(float*)malloc((size_t)maxsize/2);
  y2=(float*)malloc((size_t)maxsize/2);

  if ((y1==0)||(y2==0)) {
         ALRT("Out of memory");
      return (0);
  }

/*** Direct transform ***/

  if (nlevels>0) {

/* Computes image mean */
    xm=0.;
    kmax=(long)M*N;;
    for (kk=0;kk<kmax;kk++) xm=xm+ix[kk];
      xm=xm/(float)kmax;
      *mean=xm;

        for (k=1;k<=nlevels;k++) {
                    M2=M1/2;
                    N2=N1/2;
                  
                    pvv=wave;
                    pvw=wave+M2;
                    pwv=wave+(long)M*N2;
                    pww=pwv+M2;
                  
                  
                  /* Processing along the rows */
                            p1=ix;
                            if (k>1) p1=pvv;
                            p2=pvv;
                            nx=M1;
                            for (j=0;j<N1;j++,p1=p1+M,p2=p2+M) {
                                p11=p1;
                                if (k==1) for (i=0;i<nx;i++,p11++) x[i]=*p11-xm;
                                        else for (i=0;i<nx;i++,p11++) x[i]=*p11;
                  
                                    if (mirror==0) uqmfdownp(x,nx,y1,y2,&ny,v,nv,w,nw);
                                    else uqmfdown(x,nx,y1,y2,&ny,v,nv,w,nw);
                  
                                  p22=p2;
                                  for (i=0;i<ny;i++,p22++) *p22=y1[i];
                                  for (i=0;i<ny;i++,p22++) *p22=y2[i];
                              }
                  /* Processing along the columns */
                            p1=pvv;
                            p2=pvv;
                            p3=pwv;
                            nx=N1;
                            for (i=0;i<M2;i++,p1++,p2++,p3++) {
                                p11=p1;
                                  for (j=0;j<nx;j++,p11=p11+M) x[j]=*p11;
                  
                                  if (mirror==0) uqmfdownp(x,nx,y1,y2,&ny,v,nv,w,nw);
                                  else uqmfdown(x,nx,y1,y2,&ny,v,nv,w,nw);
                  
                                  p22=p2;p33=p3;
                                  if (k==nlevels)
                                        for (j=0;j<ny;j++,p22=p22+M,p33=p33+M) {
                                              *p22=y1[j]+xm;
                                              *p33=y2[j];
                                        }
                                  else
                                        for (j=0;j<ny;j++,p22=p22+M,p33=p33+M) {
                                              *p22=y1[j];
                                              *p33=y2[j];
                                        }
                              }
                  
                            p1=pvw;
                            p2=pvw;
                            p3=pww;
                            nx=N1;
                            for (i=0;i<M2;i++,p1++,p2++,p3++) {
                                p11=p1;
                                  for (j=0;j<nx;j++,p11=p11+M) x[j]=*p11;
                  
                                  if (mirror==0) uqmfdownp(x,nx,y1,y2,&ny,v,nv,w,nw);
                                  else uqmfdown(x,nx,y1,y2,&ny,v,nv,w,nw);
                  
                                  p22=p2;p33=p3;
                                  for (j=0;j<ny;j++,p22=p22+M,p33=p33+M) {
                                        *p22=y1[j];
                                        *p33=y2[j];
                                  }
                              }
                              M1=M2;
                              N1=N2;
        }

  } else {

/*** Indirect transform ***/

  M2=M/reduc;
  N2=N/reduc;
  xm=*mean;

      for (k=levels;k>0;k--) {
      
            M1=M2*2;
            N1=N2*2;

              pvv=wave;
              pvw=wave+M2;
              pwv=wave+(long)M*N2;
              pww=pwv+M2;
      
/* Compute the mean and substract it from lower quadrant */
#if 1
      if (k==levels) {
                  xm=0.;
                  p2=pvv;
                for (i=0;i<M2;i++,p2++) {
                      p22=p2;
                      for (j=0;j<N2;j++,p22=p22+M) xm=xm+*p22;
                  }
                  xm=xm/(float)(N2*(long)M2);
                  p2=pvv;
                for (i=0;i<M2;i++,p2++) {
                      p22=p2;
                      for (j=0;j<N2;j++,p22=p22+M) *p22=*p22-xm;
                  }
      }
#endif      
      /* Processing along the columns */
                p1=ix;
                p2=pvv;
                if (k!=levels) p2=ix;
                p3=pwv;
                ny=N2;
                for (i=0;i<M2;i++,p1++,p2++,p3++) {
                      p22=p2;p33=p3;
                      for (j=0;j<ny;j++,p22=p22+M,p33=p33+M) {
                            y1[j]=*p22;
                            y2[j]=*p33;
                      }
      
                    if (mirror==0) uqmfupp(y1,y2,ny,x,&nx,v,nv,w,nw);
                    else  uqmfup(y1,y2,ny,x,&nx,v,nv,w,nw);
      
                    p11=p1;
                      for (j=0;j<nx;j++,p11=p11+M) *p11=x[j];
      
                  }
      
                p1=ix+M2;
                p2=pvw;
                p3=pww;
                ny=N2;
                for (i=0;i<M2;i++,p1++,p2++,p3++) {
                      p22=p2;p33=p3;
                      for (j=0;j<ny;j++,p22=p22+M,p33=p33+M) {
                            y1[j]=*p22;
                            y2[j]=*p33;
                      }
      
                    if (mirror==0) uqmfupp(y1,y2,ny,x,&nx,v,nv,w,nw);
                    else uqmfup(y1,y2,ny,x,&nx,v,nv,w,nw);
      
                    p11=p1;
                      for (j=0;j<nx;j++,p11=p11+M) *p11=x[j];
      
                  }
      /* Processing along the rows */
                p1=ix;
                p2=ix;
                ny=M2;
                for (j=0;j<N1;j++,p1=p1+M,p2=p2+M) {
                      p22=p2;
                      for (i=0;i<ny;i++,p22++) y1[i]=*p22;
                      for (i=0;i<ny;i++,p22++) y2[i]=*p22;
      
                    if (mirror==0) uqmfupp(y1,y2,ny,x,&nx,v,nv,w,nw);
                    else uqmfup(y1,y2,ny,x,&nx,v,nv,w,nw);
            
                        p11=p1;
                        if (k!=1) for (i=0;i<nx;i++,p11++) *p11=x[i];
                           else for (i=0;i<nx;i++,p11++) *p11=x[i]+xm;
                  }
      
                  M2=M1;
                  N2=N1;
    }
      
      
  }      /* Indirect transform */

  free(x);
  free(y1);
  free(y2);
 return 1;
}
/*************************************************************
 *  FUNCTION : wavelet_prepare
 *  PURPOSE : Prepares data structures for representation of
 *  the wavelet transform in array storage.
 *
 *  Michael Unser / BEIP      OCT-93
 *
 *************************************************************/
void wavelet_prepare(int nx,int ny,int nz,int nlevels
                              ,int *nxw,int *nyw,int *nzw,long nxa[],long nya[],long nza[],int *nimage);
void wavelet_prepare(int nx,int ny,int nz,int nlevels
                              ,int *nxw,int *nyw,int *nzw,long nxa[],long nya[],long nza[],int *nimage)
{ int reduc,M,N,NZ,i;

       reduc=MAX(1,pow(2.,(double)nlevels));
       *nimage=1+3*nlevels;
       M=MAX(1,(nx/reduc)); N=MAX(1,(ny/reduc));  NZ=MAX(1,(nz/reduc));
       *nxw=nx > 1 ? M*reduc : 1;
       *nyw=ny > 1 ? N*reduc : 1;
       *nzw=nz > 1 ? N*reduc : 1;
       nxa[0]=M;
       nya[0]=N;
       nza[0]=NZ;
       for (i=0;i<nlevels;i++) {
              nxa[3*i+1]=M;
              nxa[3*i+2]=M;
              nxa[3*i+3]=M;
              nya[3*i+1]=N;
              nya[3*i+2]=N;
              nya[3*i+3]=N;
              nza[3*i+1]=NZ;
              nza[3*i+2]=NZ;
              nza[3*i+3]=NZ;
            if (nx>1) M=M*2;
            if (ny>1) N=N*2;
            if (nz>1) NZ=NZ*2;
      }
  return;
}
/*************************************************************
 *  FUNCTION : wav2array
 *  PURPOSE : Copies the computed wavelet image into an array of
 *  wavelet components.
 *
 *  Michael Unser / BEIP      OCT-93
 *
 *************************************************************/
void wav2array(float *wave, int nx, int ny, float *pa[], long nxa[], long
nya[], int nlevels);
void wav2array(float *wave, int nx, int ny, float *pa[], long nxa[], long
nya[], int nlevels)
{ int nimage,reduc,M,N,i,k;
  float *pvw,*pwv,*pww;

       reduc=MAX(1,pow(2.,(double)nlevels));
       nimage=1+3*nlevels;
       M=MAX(1,(nx/reduc)); N=MAX(1,(ny/reduc));
       if ((M!=nxa[0]) || (N!=nya[0])) {
         ALRT("wav2array : incorrect array dimension -> exit");
         return;
       }
       copyr(wave,nx,N,pa[0],nxa[0],nya[0]);
       for (i=0;i<nlevels;i++) {
              pvw=wave+M;
             pwv=wave+(long)nx*N;
             pww=pwv+M;
             k=3*i;
             copyr(pvw,nx,N,pa[k+1],nxa[k+1],nya[k+1]);
             copyr(pwv,nx,N,pa[k+2],nxa[k+2],nya[k+2]);
             copyr(pww,nx,N,pa[k+3],nxa[k+3],nya[k+3]);
             M=M*2;
             N=N*2;
       }
  return;
}
/*************************************************************
 *  FUNCTION : array2wave
 *  PURPOSE : Copies the array of wavelet components into a
 *  single float image.
 *
 *  Michael Unser / BEIP      OCT-93
 *
 *************************************************************/
void array2wav(float *wave, int nx, int ny, float *pa[], long nxa[], long
nya[], int nlevels);
void array2wav(float *wave, int nx, int ny, float *pa[], long nxa[], long
nya[], int nlevels)
{ int nimage,reduc,M,N,i,k;
  float *pvw,*pwv,*pww;

       reduc=MAX(1,pow(2.,(double)nlevels));
       nimage=1+3*nlevels;
       M=MAX(1,(nx/reduc)); N=MAX(1,(ny/reduc));
       if ((M!=nxa[0]) || (N!=nya[0])) {
         ALRT("array2wave : incorrect array dimension -> exit");
         return;
       }
       copyr(pa[0],nxa[0],nya[0],wave,nx,N);
       for (i=0;i<nlevels;i++) {
              pvw=wave+M;
             pwv=wave+(long)nx*N;
             pww=pwv+M;
             k=3*i;
             copyr(pa[k+1],nxa[k+1],nya[k+1],pvw,nx,N);
             copyr(pa[k+2],nxa[k+2],nya[k+2],pwv,nx,N);
             copyr(pa[k+3],nxa[k+3],nya[k+3],pww,nx,N);
             M=M*2;
             N=N*2;
       }
  return;
}
/*************************************************************
 *  FUNCTION : threshold_channel
 *  PURPOSE : Performs thresholding in wavelet domain for signal
 *  denoising.
 *
 *      PARAMETERS : t = threshold, level0 = finer level of resolution
 *
 *  Michael Unser / BEIP      AUG-94
 *
 *************************************************************/
      void threshold_channel(float *p, int nx, int nx0, int ny, float t, long
*icount);

void wave_threshold(float *wave, int nx, int ny, int nlevels,  int level0,
float t0, long *icount);
void wave_threshold(float *wave, int nx, int ny, int nlevels,  int level0,
float t0, long *icount)
{ int nimage,reduc,M,N,i,k;
  float *pvw,*pwv,*pww;
  float t;

       *icount=0;
       reduc=MAX(1,pow(2.,(double)nlevels));
       nimage=1+3*nlevels;
       M=MAX(1,(nx/reduc)); N=MAX(1,(ny/reduc));
       t=t0/(float)reduc;

       threshold_channel(wave,M,nx,N,t,icount);
       for (i=0;i<nlevels;i++) {
              pvw=wave+M;
             pwv=wave+(long)nx*N;
             pww=pwv+M;
             k=3*i;
             if (i>(nlevels-level0)) t=1000.E10;            /* all channels finer than level0
are set to zero */

             threshold_channel(pvw,M,nx,N,t,icount);
             threshold_channel(pwv,M,nx,N,t,icount);
             threshold_channel(pww,M,nx,N,t,icount);
             M=M*2;
             N=N*2;
             t=t*2;
       }
}
void threshold_channel(float *p, int nx, int nx0, int ny, float t, long
*icount);
void threshold_channel(float *p, int nx, int nx0, int ny, float t, long *icount)
{ float *p0;
  int i,j;
  for (j=0;j<ny;j++) {
    p0=p+(long)nx0*(long)j;
    for (i=0;i<nx;i++) {
      if (fabs(p0[i])<=t) p0[i]=0.;
      else *icount=*icount+1;

    }
  }

}
/*************************************************************
 *  FUNCTION : uqmfdown
 *  PURPOSE : Two channel filtering and down-sampling for 1D
 *  real signals. The filters are symmetrical and the signal
 *      is extended using mirror extrapolation.
 *
 *      lowpass filter : v[nv], highpass filter : w[nw]
 *
 *      x[nx] -> y1[nx/2] & y2[nx/2]
 *
 *  Michael Unser / BEIP      OCT-92      (tested)
 *
 *************************************************************/
void uqmfdown(float x[],int nx,float y1[],float y2[],int *ny,float v[],int
nv,float w[],int nw);
void uqmfdown(float x[],int nx,float y1[],float y2[],int *ny,float v[],int
nv,float w[],int nw)
{        float pix;
        int i,j,k,nyy,n,j1,j2,kn;

        nyy=nx/2;
        n=nyy*2;                                    /* signal length */
        kn=2*n-2;                                    /* Global period */

  if ((nv<=1)||(nw<=1))                        /* Haar transform */
    for (i=0;i<nyy;i++) {
      j=2*i;
      y1[i]=(+x[j]+x[j+1])/2.;
      y2[i]=(-x[j]+x[j+1])/2.;
      }
  else
    for (i=0;i<nyy;i++) {
      j=i*2;
      pix=x[j]*v[0];
      for (k=1;k<nv;k++) {
            j1=j-k;
            j2=j+k;
            
            if (j1<0) {
              j1=(-j1) % kn;
              if (j1>n-1) j1=kn-j1;
            }
            if (j2>n-1) {
              j2=j2 % kn;
              if (j2>n-1) j2=kn-j2;
            }

            pix=pix+v[k]*(x[j1]+x[j2]);
      }
      y1[i]=pix;

      j=j+1;
      pix=x[j]*w[0];
      for (k=1;k<nw;k++) {
            j1=j-k;
            j2=j+k;

            if (j1<0) {
              j1=(-j1) % kn;
              if (j1>n-1) j1=kn-j1;
            }
            if (j2>n-1) {
              j2=j2 % kn;
              if (j2>n-1) j2=kn-j2;
            }

            pix=pix+w[k]*(x[j1]+x[j2]);
      }
      y2[i]=pix;
    }
      *ny=nyy;
}
/*************************************************************
 *  FUNCTION : uqmfup
 *  PURPOSE : Two channel up-sampling and filtering for 1D
 *  real signals. The filters are symmetrical and the
 *  boundary conditions are compatible with uqmfdown().
 *
 *      lowpass filter : v[nv], highpass filter : w[nw]
 *
 *      y1[ny] & y2[ny] -> x[2*ny]
 *
 *  Michael Unser / BEIP      OCT-92      (tested)
 *                          APR-97      (corrected boundary conditions)
 *
 *************************************************************/

void uqmfup(float y1[],float y2[],int ny,float x[],int *nx,float v[],int
nv,float w[],int nw);
void uqmfup(float y1[],float y2[],int ny,float x[],int *nx,float v[],int
nv,float w[],int nw)

{        float pix1,pix2;
        int i,j,k,kk,i1,i2,k01,k02,n,kn;

        *nx=ny*2;
        k01=(nv/2)*2-1;
        k02=(nw/2)*2-1;
        n=ny;
        kn=2*n-1;

      if ((nv<=1)||(nw<=1))                        /* Haar transform */
          for (i=0;i<ny;i++) {
            j=2*i;
            x[j]=y1[i]-y2[i];
            x[j+1]=y1[i]+y2[i];
            }
      else
      for (i=0;i<ny;i++)      {
            j=2*i;
            pix1=v[0]*y1[i];
            for (k=2;k<nv;k=k+2) {
              i1=i-(k/2);
              i2=i+(k/2);
                  if (i1<0) {                        /* standard boundary */
                    i1=(-i1) % kn;
                    if (i1>n-1) i1=kn-i1;
                  }
                  if (i2>n-1) {                  /* non-standard boundary */
                    i2=(i2) % kn;
                    if (i2>n-1) i2=kn-i2;
                  }
/* NOTE : This is the correct boundary condition assuming that the original
   image has an even size (2n) and the decimation was performed according to
   the sequence 1,3,5,I,2n-1.      */
              pix1=pix1+v[k]*(y1[i1]+y1[i2]);
            }
      
            pix2=0.;
            for (k=-k02;k<nw;k=k+2) {
                  kk=abs(k);
                  i1=i+(k-1)/2;
              if (i1<0) {                        /* non-standard boundary */
                    i1=(-i1-1) % kn;
                    if (i1>n-1) i1=kn-1-i1;
              }
                  if (i1>n-1) {                  /* standard boundary */
                    i1=i1 % kn;
                    if (i1>n-1) i1=kn-1-i1;
                  }
/* NOTE : This is the correct boundary condition assuming that the original
   image has an even size (2n) and the decimation was performed according to
   the sequence 2,4,6,I,2n.      */
              pix2=pix2+w[kk]*y2[i1];
            }
/*            x[j]=(pix1+pix2)*2;      */
            x[j]=(pix1+pix2);

            j=j+1;
            pix1=0.;
            for (k=-k01;k<nv;k=k+2) {
                      kk=abs(k);
                        i1=i+(k+1)/2;
                        if (i1<0) {                        /* standard boundary */
                          i1=(-i1) % kn;
                          if (i1>n-1) i1=kn-i1;
                        }
                        if (i1>n-1) {                  /* non-standard boundary */
                          i1=(i1) % kn;
                          if (i1>n-1) i1=kn-i1;
                        }
              pix1=pix1+v[kk]*y1[i1];
            }
            pix2=w[0]*y2[i];
            for (k=2;k<nw;k=k+2) {
              i1=i-(k/2);
              i2=i+(k/2);
              if (i1<0) {                        /* non-standard boundary */
                    i1=(-i1-1) % kn;
                    if (i1>n-1) i1=kn-1-i1;
              }
                  if (i2>n-1) {                  /* standard boundary */
                    i2=i2 % kn;
                    if (i2>n-1) i2=kn-1-i2;
                  }
              pix2=pix2+w[k]*(y2[i1]+y2[i2]);
            }
/*            x[j]=(pix1+pix2)*2;      */
              x[j]=(pix1+pix2);
      }
}
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 *  PERIODIC BOUNDARY CONDITIONS :
 * +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/*************************************************************
 *  FUNCTION : uqmfdownp
 *  PURPOSE : Two channel filtering and down-sampling for 1D
 *  real signals. The filters are symmetrical and the signal
 *      is extended using periodic boundary conditions.
 *
 *      lowpass filter : v[nv], highpass filter : w[nw]
 *
 *      x[nx] -> y1[nx/2] & y2[nx/2]
 *
 *  Michael Unser / BEIP      OCT-92      (tested)
 *
 *************************************************************/
void uqmfdownp(float x[],int nx,float y1[],float y2[],int *ny,float v[],int
nv,float w[],int nw);
void uqmfdownp(float x[],int nx,float y1[],float y2[],int *ny,float v[],int nv,float w[],int nw)
{        float pix;
        int i,j,k,nyy,n,j1,j2,kn,knn;

        nyy=nx/2;
        n=nyy*2;                                    /* signal length */
        kn=n;                                        /* Global period */
        knn=20*n;

  if ((nv<=1)||(nw<=1))                        /* Haar transform */
    for (i=0;i<nyy;i++) {
      j=2*i;
      y1[i]=(+x[j]+x[j+1])/2.;
      y2[i]=(-x[j]+x[j+1])/2.;
      }
  else
    for (i=0;i<nyy;i++) {
      j=i*2;
      pix=x[j]*v[0];
      for (k=1;k<nv;k++) {
            j1=j-k;
            j2=j+k;
            
            if (j1<0) j1=(knn+j1) % kn;
            if (j2>n-1) j2=j2 % kn;

            pix=pix+v[k]*(x[j1]+x[j2]);
      }
      y1[i]=pix;

      j=j+1;
      pix=x[j]*w[0];
      for (k=1;k<nw;k++) {
            j1=j-k;
            j2=j+k;

            if (j1<0) j1=(knn+j1) % kn;
            if (j2>n-1) j2=j2 % kn;

            pix=pix+w[k]*(x[j1]+x[j2]);
      }
      y2[i]=pix;
    }
      *ny=nyy;
}
/*************************************************************
 *  FUNCTION : uqmfupp
 *  PURPOSE : Two channel up-sampling and filtering for 1D
 *  real signals. The filters are symmetrical and the
 *  boundary conditions are compatible with uqmfdownp().
 *
 *      lowpass filter : v[nv], highpass filter : w[nw]
 *
 *      y1[ny] & y2[ny] -> x[2*ny]
 *
 *  Michael Unser / BEIP      OCT-92      (tested)
 *
 *************************************************************/

void uqmfupp(float y1[],float y2[],int ny,float x[],int *nx,float v[],int
nv,float w[],int nw);
void uqmfupp(float y1[],float y2[],int ny,float x[],int *nx,float v[],int
nv,float w[],int nw)

{        float pix1,pix2;
        int i,j,k,kk,i1,i2,k01,k02,n,kn,knn;

        *nx=ny*2;
        k01=(nv/2)*2-1;
        k02=(nw/2)*2-1;
        n=ny;
        kn=n;
        knn=20*n;

      if ((nv<=1)||(nw<=1))                        /* Haar transform */
          for (i=0;i<ny;i++) {
            j=2*i;
            x[j]=y1[i]-y2[i];
            x[j+1]=y1[i]+y2[i];
            }
      else
      for (i=0;i<ny;i++)      {
            j=2*i;
            pix1=v[0]*y1[i];
            for (k=2;k<nv;k=k+2) {
              i1=i-(k/2);
              i2=i+(k/2);

                  if (i1<0) i1=(knn+i1) % kn;
                  if (i2>n-1) i2=i2 % kn;

              pix1=pix1+v[k]*(y1[i1]+y1[i2]);
            }
      
            pix2=0.;
            for (k=-k02;k<nw;k=k+2) {
                  kk=abs(k);
                  i1=i+(k-1)/2;
              if (i1<0) i1=(knn+i1) % kn;
                  if (i1>n-1) i1=i1 % kn;

              pix2=pix2+w[kk]*y2[i1];
            }
            x[j]=(pix1+pix2);

            j=j+1;
            pix1=0.;
            for (k=-k01;k<nv;k=k+2) {
                      kk=abs(k);
                        i1=i+(k+1)/2;
                        if (i1<0) i1=(knn+i1) % kn;
                        if (i1>n-1) i1=i1 % kn;

              pix1=pix1+v[kk]*y1[i1];
            }
            pix2=w[0]*y2[i];
            for (k=2;k<nw;k=k+2) {
              i1=i-(k/2);
              i2=i+(k/2);
              if (i1<0) i1=(knn+i1) % kn;
                  if (i2>n-1) i2=i2 % kn;

              pix2=pix2+w[k]*(y2[i1]+y2[i2]);
            }
              x[j]=(pix1+pix2);
      }
}




ASKER CERTIFIED SOLUTION
Avatar of atlan2
atlan2
Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
Avatar of AdamB_Open
AdamB_Open

ASKER

Comment accepted as answer
Avatar of AdamB_Open
AdamB_Open

ASKER

Thanks,
   this is exactly what I'm after.  I haven't tested it yet but it looks like what I need.
Adam.