[Last Call] Learn about multicloud storage options and how to improve your company's cloud strategy. Register Now

x
?
Solved

Wavelets

Posted on 2000-04-13
5
Medium Priority
?
212 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 1200 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

Technology Partners: We Want Your Opinion!

We value your feedback.

Take our survey and automatically be enter to win anyone of the following:
Yeti Cooler, Amazon eGift Card, and Movie eGift Card!

Question has a verified solution.

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

A short article about problems I had with the new location API and permissions in Marshmallow
We live in a world of interfaces like the one in the title picture. VBA also allows to use interfaces which offers a lot of possibilities. This article describes how to use interfaces in VBA and how to work around their bugs.
Simple Linear Regression
Starting up a Project

650 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