Solved

Wavelets

Posted on 2000-04-13
5
203 Views
Last Modified: 2010-04-16
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.
0
Comment
Question by:AdamB_Open
  • 2
  • 2
5 Comments
 
LVL 22

Expert Comment

by:cookre
Comment Utility
If you don't get any direct responses, you may want to post some of the pseudo code for folks here to work on.
0
 
LVL 1

Expert Comment

by:atlan2
Comment Utility
"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);
      }
}




0
 
LVL 1

Accepted Solution

by:
atlan2 earned 300 total points
Comment Utility
"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);
      }
}




0
 

Author Comment

by:AdamB_Open
Comment Utility
Comment accepted as answer
0
 

Author Comment

by:AdamB_Open
Comment Utility
Thanks,
   this is exactly what I'm after.  I haven't tested it yet but it looks like what I need.
Adam.
0

Featured Post

How to improve team productivity

Quip adds documents, spreadsheets, and tasklists to your Slack experience
- Elevate ideas to Quip docs
- Share Quip docs in Slack
- Get notified of changes to your docs
- Available on iOS/Android/Desktop/Web
- Online/Offline

Join & Write a Comment

Suggested Solutions

Title # Comments Views Activity
 shows up in Outlook, not OWA or on phone 3 57
powerN  challenge 3 46
Specific format 21 139
Device same like our heart 12 42
Here we come across an interesting topic of coding guidelines while designing automation test scripts. The scope of this article will not be limited to QTP but to an overall extent of using VB Scripting for automation projects. Introduction Now…
Although it can be difficult to imagine, someday your child will have a career of his or her own. He or she will likely start a family, buy a home and start having their own children. So, while being a kid is still extremely important, it’s also …
In this fourth video of the Xpdf series, we discuss and demonstrate the PDFinfo utility, which retrieves the contents of a PDF's Info Dictionary, as well as some other information, including the page count. We show how to isolate the page count in a…
In this seventh video of the Xpdf series, we discuss and demonstrate the PDFfonts utility, which lists all the fonts used in a PDF file. It does this via a command line interface, making it suitable for use in programs, scripts, batch files — any pl…

728 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question

Need Help in Real-Time?

Connect with top rated Experts

9 Experts available now in Live!

Get 1:1 Help Now