# Wavelets

I need an algorith to do a forward wavelet transformation of 2-D data using Gabor wavelets.  The algorithm will be implemented in a specialised language, but I can read VB,C/C++.  I can find plenty of mathematical explanations and pseudocode for wavelets on the net, but I have difficulty converting this into code.
###### Who is Participating?

x

Commented:
"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

Commented:
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

Commented:
"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 Commented:
0

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