Solved

Wavelets

Posted on 2000-04-13
5
210 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
[X]
Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

  • Help others & share knowledge
  • Earn cash & points
  • Learn & ask questions
  • 2
  • 2
5 Comments
 
LVL 22

Expert Comment

by:cookre
ID: 2717465
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
ID: 2844792
"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
ID: 2844817
"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
ID: 2847694
Comment accepted as answer
0
 

Author Comment

by:AdamB_Open
ID: 2847695
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

[Webinar] Learn How Hackers Steal Your Credentials

Do You Know How Hackers Steal Your Credentials? Join us and Skyport Systems to learn how hackers steal your credentials and why Active Directory must be secure to stop them. Thursday, July 13, 2017 10:00 A.M. PDT

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

This article will show, step by step, how to integrate R code into a R Sweave document
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 …
An introduction to basic programming syntax in Java by creating a simple program. Viewers can follow the tutorial as they create their first class in Java. Definitions and explanations about each element are given to help prepare viewers for future …
Introduction to Processes

691 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