telafft.c

用于Pattern Recongnition的一个图像处理程序－－Trace Transform。

    nf=ifac[1];
    na=0;
    l1=1;
    iw=0;
    for(k1=1; k1<=nf; k1++)
    {
	ip=ifac[k1+1];
	l2=ip*l1;
	ido=n / l2;
	idl1=ido*l1;
	if(ip==4)
	{
	    ix2=iw+ido;
	    ix3=ix2+ido;
	    if(na==0)
		radb4(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
	    else
		radb4(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
	    na=1-na;
	}
	else if(ip==2)
	{
	    if(na==0)
		radb2(ido, l1, c, ch, &wa[iw]);
	    else
		radb2(ido, l1, ch, c, &wa[iw]);
	    na=1-na;
	}
	else if(ip==3)
	{
	    ix2=iw+ido;
	    if(na==0)
		radb3(ido, l1, c, ch, &wa[iw], &wa[ix2]);
	    else
		radb3(ido, l1, ch, c, &wa[iw], &wa[ix2]);
	    na=1-na;
	}
	else if(ip==5)
	{
	    ix2=iw+ido;
	    ix3=ix2+ido;
	    ix4=ix3+ido;
	    if(na==0)
		radb5(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
	    else
		radb5(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
	    na=1-na;
	}
	else
	{
	    if(na==0)
		radbg(ido, ip, l1, idl1, c, c, c, ch, ch, &wa[iw]);
	    else
		radbg(ido, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw]);
	    if(ido==1)
		na=1-na;
	}
	l1=l2;
	iw+=(ip-1)*ido;
    }
    if(na==0)
	return;
    for(i=0; i<n; i++)
	c[i]=ch[i];
}				/*rfftb1*/


void 
rfftf(int n, double r[], double wsave[])
{
    if(n==1)
	return;
    rfftf1(n, r, wsave, wsave+n,(int*)(wsave+2*n));
}				/*rfftf*/


void 
rfftb(int n, double r[], double wsave[])
{
    if(n==1)
	return;
    rfftb1(n, r, wsave, wsave+n,(int*)(wsave+2*n));
}				/*rfftb*/


static void 
rffti1(int n, double wa[], int ifac[])
{
    static const int ntryh[4]=
    {4, 2, 3, 5};
    static const double twopi=6.28318530717959;
    double  argh;
    int     ntry, i, j;
    double  argld;
    int     k1, l1, l2, ib;
    double  fi;
    int     ld, ii, nf, ip, nl, is, nq, nr;
    double  arg;
    int     ido, ipm;
    int     nfm1;

    nl=n;
    nf=0;
    j=0;
  startloop:
   ++j;
    if(j<=4)
	ntry=ntryh[j-1];
    else
	ntry+=2;
    do
    {
	nq=nl / ntry;
	nr=nl-ntry*nq;
	if(nr !=0)
	    goto startloop;
	++nf;
	ifac[nf+1]=ntry;
	nl=nq;
	if(ntry==2 && nf !=1)
	{
	    for(i=2; i<=nf; i++)
	    {
		ib=nf-i+2;
		ifac[ib+1]=ifac[ib];
	    }
	    ifac[2]=2;
	}
    }
    while(nl !=1);
    ifac[0]=n;
    ifac[1]=nf;
    argh=twopi /(double)(n);
    is=0;
    nfm1=nf-1;
    l1=1;
    if(nfm1==0)
	return;
    for(k1=1; k1<=nfm1; k1++)
    {
	ip=ifac[k1+1];
	ld=0;
	l2=l1*ip;
	ido=n / l2;
	ipm=ip-1;
	for(j=1; j<=ipm;++j)
	{
	    ld+=l1;
	    i=is;
	    argld=(double)ld*argh;

	    fi=0;
	    for(ii=3; ii<=ido; ii+=2)
	    {
		i+=2;
		fi+=1;
		arg=fi*argld;
		wa[i-2]=cos(arg);
		wa[i-1]=sin(arg);
	    }
	    is+=ido;
	}
	l1=l2;
    }
}				/*rffti1*/


void 
rffti(int n, double wsave[])
{
    if(n==1)
	return;
    rffti1(n, wsave+n,(int*)(wsave+2*n));
}				/*rffti*/

/*----------------------------------------------------------------------
   cosqf1, cosqb1, cosqf, cosqb, cosqi. Quarter-wave cos-FFTs.
  ----------------------------------------------------------------------*/

static void 
cosqf1(int n, double x[], double w[], double xh[])
{
    int     modn, i, k;
    int     kc, np2, ns2;
    double  xim1;

    ns2=(n+1)/ 2;
    np2=n+2;
    for(k=1; k<ns2; k++)
    {
	kc=n-k;
	xh[k]=x[k]+x[kc];
	xh[kc]=x[k]-x[kc];
    }
    modn=n%2;
    if(modn==0)
	xh[ns2]=x[ns2]+x[ns2];
    for(k=1; k<ns2; k++)
    {
	kc=n-k;
	x[k]=w[k-1]*xh[kc]+w[kc-1]*xh[k];
	x[kc]=w[k-1]*xh[k]-w[kc-1]*xh[kc];
    }
    if(modn==0)
	x[ns2]=w[ns2-1]*xh[ns2];
    rfftf(n, x, xh);
    for(i=2; i<n; i+=2)
    {
	xim1=x[i-1]-x[i];
	x[i]=x[i-1]+x[i];
	x[i-1]=xim1;
    }
}				/*cosqf1*/


static void 
cosqb1(int n, double x[], double w[], double xh[])
{
    int     modn, i, k;
    int     kc, ns2;
    double  xim1;

    ns2=(n+1)/ 2;
    for(i=2; i<n; i+=2)
    {
	xim1=x[i-1]+x[i];
	x[i]-=x[i-1];
	x[i-1]=xim1;
    }
    x[0]+=x[0];
    modn=n%2;
    if(modn==0)
	x[n-1]+=x[n-1];
    rfftb(n, x, xh);
    for(k=1; k<ns2; k++)
    {
	kc=n-k;
	xh[k]=w[k-1]*x[kc]+w[kc-1]*x[k];
	xh[kc]=w[k-1]*x[k]-w[kc-1]*x[kc];
    }
    if(modn==0)
	x[ns2]=w[ns2-1]*(x[ns2]+x[ns2]);
    for(k=1; k<ns2; k++)
    {
	kc=n-k;
	x[k]=xh[k]+xh[kc];
	x[kc]=xh[k]-xh[kc];
    }
    x[0]+=x[0];
}				/*cosqb1*/


void 
cosqf(int n, double x[], double wsave[])
{
    static const double sqrt2=1.4142135623731;
    double  tsqx;

    if(n<2)
    {
	return;
    }
    else if(n==2)
    {
	tsqx=sqrt2*x[1];
	x[1]=x[0]-tsqx;
	x[0]+=tsqx;
    }
    else
    {
	cosqf1(n, x, wsave, &wsave[n]);
    }
}				/*cosqf*/


void 
cosqb(int n, double x[], double wsave[])
{
    static const double tsqrt2=2.82842712474619;
    double  x1;

    if(n<2)
    {
	x[0]*=4;
    }
    else if(n==2)
    {
	x1=4*(x[0]+x[1]);
	x[1]=tsqrt2*(x[0]-x[1]);
	x[0]=x1;
    }
    else
	cosqb1(n, x, wsave, &wsave[n]);
}				/*cosqb*/


void 
cosqi(int n, double wsave[])
{
    static const double pih=1.57079632679491;
    int     k;
    double  dt;

    dt=pih / n;
    for(k=0; k<n; k++)
	wsave[k]=cos((k+1)*dt);
    rffti(n, &wsave[n]);
}				/*cosqi*/

/*----------------------------------------------------------------------
   cost, costi. Full-wave cos-FFTs.
  ----------------------------------------------------------------------*/

void 
cost(int n, double x[], double wsave[])
{
    int     modn, i, k;
    double  c1, t1, t2;
    int     kc;
    double  xi;
    int     nm1;
    double  x1h;
    int     ns2;
    double  tx2, x1p3, xim2;

    nm1=n-1;
    ns2=n / 2;
    if(n-2<0)
	return;
    else if(n==2)
    {
	x1h=x[0]+x[1];
	x[1]=x[0]-x[1];
	x[0]=x1h;
    }
    else if(n==3)
    {
	x1p3=x[0]+x[2];
	tx2=x[1]+x[1];
	x[1]=x[0]-x[2];
	x[0]=x1p3+tx2;
	x[2]=x1p3-tx2;
    }
    else
    {
	c1=x[0]-x[n-1];
	x[0]+=x[n-1];
	for(k=1; k<ns2; k++)
	{
	    kc=nm1-k;
	    t1=x[k]+x[kc];
	    t2=x[k]-x[kc];
	    c1+=wsave[kc]*t2;
	    t2=wsave[k]*t2;
	    x[k]=t1-t2;
	    x[kc]=t1+t2;
	}
	modn=n%2;
	if(modn !=0)
	    x[ns2]+=x[ns2];
	rfftf(nm1, x, &wsave[n]);
	xim2=x[1];
	x[1]=c1;
	for(i=3; i<n; i+=2)
	{
	    xi=x[i];
	    x[i]=x[i-2]-x[i-1];
	    x[i-1]=xim2;
	    xim2=xi;
	}
	if(modn !=0)
	    x[n-1]=xim2;
    }
}				/*cost*/


void 
costi(int n, double wsave[])
{
    static const double pi=3.14159265358979;
    int     k, kc, ns2;
    double  dt;

    if(n<=3)
	return;
    ns2=n / 2;
    dt=pi /(n-1);
    for(k=1; k<ns2; k++)
    {
	kc=n-k-1;
	wsave[k]=2*sin(k*dt);
	wsave[kc]=2*cos(k*dt);
    }
    rffti(n-1, &wsave[n]);
}				/*costi*/

/*----------------------------------------------------------------------
   sinqf, sinqb, sinqi. Quarter-wave sin-FFTs. These call cosqf etc.
  ----------------------------------------------------------------------*/

void 
sinqf(int n, double x[], double wsave[])
{
    int     k;
    double  xhold;
    int     kc, ns2;

    if(n==1)
	return;
    ns2=n / 2;
    for(k=0; k<ns2; k++)
    {
	kc=n-k-1;
	xhold=x[k];
	x[k]=x[kc];
	x[kc]=xhold;
    }
    cosqf(n, x, wsave);
    for(k=1; k<n; k+=2)
	x[k]=-x[k];
}				/*sinqf*/


void 
sinqb(int n, double x[], double wsave[])
{
    int     k;
    double  xhold;
    int     kc, ns2;

    if(n<=1)
    {
	x[0]*=4;
	return;
    }
    ns2=n / 2;
    for(k=1; k<n; k+=2)
	x[k]=-x[k];
    cosqb(n, x, wsave);
    for(k=0; k<ns2; k++)
    {
	kc=n-k-1;
	xhold=x[k];
	x[k]=x[kc];
	x[kc]=xhold;
    }
}				/*sinqb*/


void 
sinqi(int n, double wsave[])
{
    cosqi(n, wsave);
}

/*----------------------------------------------------------------------
   sint1, sint, sinti. Full-wave sin-FFTs.
  ----------------------------------------------------------------------*/

static void 
sint1(int n, double war[], const double was[], double xh[], double x[], const int ifac[])
{
    static const double sqrt3=1.73205080756888;
    int     modn, i, k;
    double  xhold, t1, t2;
    int     kc, np1, ns2;

    for(i=0; i<n; i++)
    {
	xh[i]=war[i];
	war[i]=x[i];
    }
    if(n<2)
    {
	xh[0]+=xh[0];
    }
    else if(n==2)
    {
	xhold=sqrt3*(xh[0]+xh[1]);
	xh[1]=sqrt3*(xh[0]-xh[1]);
	xh[0]=xhold;
    }
    else
    {
	np1=n+1;
	ns2=n / 2;
	x[0]=0;
	for(k=0; k<ns2; k++)
	{
	    kc=n-k-1;
	    t1=xh[k]-xh[kc];
	    t2=was[k]*(xh[k]+xh[kc]);
	    x[k+1]=t1+t2;
	    x[kc+1]=t2-t1;
	}
	modn=n%2;
	if(modn !=0)
	    x[ns2+1]=4*xh[ns2];
	rfftf1(np1, x, xh, war, ifac);
	xh[0]=0.5*x[0];
	for(i=2; i<n; i+=2)
	{
	    xh[i-1]=-x[i];
	    xh[i]=xh[i-2]+x[i-1];
	}
	if(modn==0)
	    xh[n-1]=-x[n];
    }
    for(i=0; i<n; i++)
    {
	x[i]=war[i];
	war[i]=xh[i];
    }
}				/*sint1*/


void 
sint(int n, double x[], double wsave[])
{
    int     iw1, iw2, iw3;

    iw1=n / 2;
    iw2=iw1+n+1;
    iw3=iw2+n+1;
    sint1(n, x, wsave, wsave+iw1, wsave+iw2,(int*)(wsave+iw3));
}				/*sint*/


void 
sinti(int n, double wsave[])
{
    static const double pi=3.14159265358979;
    int     k, ns2;
    double  dt;

    if(n<=1)
	return;
    ns2=n / 2;
    dt=pi /(n+1);
    for(k=0; k<ns2; k++)
	wsave[k]=2*sin((k+1)*dt);
    rffti(n+1, &wsave[ns2]);
}				/*sinti*/