Solved

# Wavelets

Posted on 2000-04-13
Medium Priority
215 Views
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
• 2
• 2

LVL 22

Expert Comment

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

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

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

ID: 2847694
0

Author Comment

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

## Featured Post

Question has a verified solution.

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

When you see single cell contains number and text, and you have to get any date out of it seems like cracking our heads.
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.
In this fifth video of the Xpdf series, we discuss and demonstrate the PDFdetach utility, which is able to list and, more importantly, extract attachments that are embedded in PDF files. It does this via a command line interface, making it suitable …
In this seventh video of the Xpdf series, we discuss and demonstrate the PDFfonts utility, which lists all the fonts used in a PDF file. It does this via a command line interface, making it suitable for use in programs, scripts, batch files — any pl…
###### Suggested Courses
Course of the Month6 days, 6 hours left to enroll