telafft.c

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

		ch[(k+jc*l1)*ido];
	}
    }
    if(ido==1)
	return;
    if(nbd>=l1)
    {
	for(j=1; j<ipph; j++)
	{
	    jc=ip-j;
	    j2=2*j;
	    for(k=0; k<l1; k++)
	    {
		for(i=2; i<ido; i+=2)
		{
		    ic=ido-i;
		    cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
		    cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
		    cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
		    cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
		}
	    }
	}
    }
    else
    {
	for(j=1; j<ipph; j++)
	{
	    jc=ip-j;
	    j2=2*j;
	    for(i=2; i<ido; i+=2)
	    {
		ic=ido-i;
		for(k=0; k<l1; k++)
		{
		    cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
		    cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
		    cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
		    cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
		}
	    }
	}
    }
}				/*radfg*/


static void 
radbg(int ido, int ip, int l1, int idl1,
       double cc[], double c1[], double c2[], double ch[], double ch2[], const double wa[])
{
    static const double twopi=6.28318530717959;
    int     idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
    double  dc2, ai1, ai2, ar1, ar2, ds2;
    int     nbd;
    double  dcp, arg, dsp, ar1h, ar2h;

    arg=twopi / ip;
    dcp=cos(arg);
    dsp=sin(arg);
    nbd=(ido-1)/ 2;
    ipph=(ip+1)/ 2;
    if(ido>=l1)
    {
	for(k=0; k<l1; k++)
	{
	    for(i=0; i<ido; i++)
	    {
		ch[i+k*ido]=cc[i+k*ip*ido];
	    }
	}
    }
    else
    {
	for(i=0; i<ido; i++)
	{
	    for(k=0; k<l1; k++)
	    {
		ch[i+k*ido]=cc[i+k*ip*ido];
	    }
	}
    }
    for(j=1; j<ipph; j++)
    {
	jc=ip-j;
	j2=2*j;
	for(k=0; k<l1; k++)
	{
	    ch[(k+j*l1)*ido]=cc[ido-1+(j2-1+k*ip)*ido]+cc[ido-1+(j2-1+k*ip)*ido];
	    ch[(k+jc*l1)*ido]=cc[(j2+k*ip)*ido]+cc[(j2+k*ip)*ido];
	}
    }

    if(ido !=1)
    {
	if(nbd>=l1)
	{
	    for(j=1; j<ipph; j++)
	    {
		jc=ip-j;
		for(k=0; k<l1; k++)
		{
		    for(i=2; i<ido; i+=2)
		    {
			ic=ido-i;
			ch[i-1+(k+j*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
			ch[i-1+(k+jc*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
			ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
			ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
		    }
		}
	    }
	}
	else
	{
	    for(j=1; j<ipph; j++)
	    {
		jc=ip-j;
		for(i=2; i<ido; i+=2)
		{
		    ic=ido-i;
		    for(k=0; k<l1; k++)
		    {
			ch[i-1+(k+j*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
			ch[i-1+(k+jc*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
			ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
			ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
		    }
		}
	    }
	}
    }

    ar1=1;
    ai1=0;
    for(l=1; l<ipph; l++)
    {
	lc=ip-l;
	ar1h=dcp*ar1-dsp*ai1;
	ai1=dcp*ai1+dsp*ar1;
	ar1=ar1h;
	for(ik=0; ik<idl1; ik++)
	{
	    c2[ik+l*idl1]=ch2[ik]+ar1*ch2[ik+idl1];
	    c2[ik+lc*idl1]=ai1*ch2[ik+(ip-1)*idl1];
	}
	dc2=ar1;
	ds2=ai1;
	ar2=ar1;
	ai2=ai1;
	for(j=2; j<ipph; j++)
	{
	    jc=ip-j;
	    ar2h=dc2*ar2-ds2*ai2;
	    ai2=dc2*ai2+ds2*ar2;
	    ar2=ar2h;
	    for(ik=0; ik<idl1; ik++)
	    {
		c2[ik+l*idl1]+=ar2*ch2[ik+j*idl1];
		c2[ik+lc*idl1]+=ai2*ch2[ik+jc*idl1];
	    }
	}
    }
    for(j=1; j<ipph; j++)
    {
	for(ik=0; ik<idl1; ik++)
	{
	    ch2[ik]+=ch2[ik+j*idl1];
	}
    }
    for(j=1; j<ipph; j++)
    {
	jc=ip-j;
	for(k=0; k<l1; k++)
	{
	    ch[(k+j*l1)*ido]=c1[(k+j*l1)*ido]-c1[(k+jc*l1)*ido];
	    ch[(k+jc*l1)*ido]=c1[(k+j*l1)*ido]+c1[(k+jc*l1)*ido];
	}
    }

    if(ido==1)
	return;
    if(nbd>=l1)
    {
	for(j=1; j<ipph; j++)
	{
	    jc=ip-j;
	    for(k=0; k<l1; k++)
	    {
		for(i=2; i<ido; i+=2)
		{
		    ch[i-1+(k+j*l1)*ido]=c1[i-1+(k+j*l1)*ido]-c1[i+(k+jc*l1)*ido];
		    ch[i-1+(k+jc*l1)*ido]=c1[i-1+(k+j*l1)*ido]+c1[i+(k+jc*l1)*ido];
		    ch[i+(k+j*l1)*ido]=c1[i+(k+j*l1)*ido]+c1[i-1+(k+jc*l1)*ido];
		    ch[i+(k+jc*l1)*ido]=c1[i+(k+j*l1)*ido]-c1[i-1+(k+jc*l1)*ido];
		}
	    }
	}
    }
    else
    {
	for(j=1; j<ipph; j++)
	{
	    jc=ip-j;
	    for(i=2; i<ido; i+=2)
	    {
		for(k=0; k<l1; k++)
		{
		    ch[i-1+(k+j*l1)*ido]=c1[i-1+(k+j*l1)*ido]-c1[i+(k+jc*l1)*ido];
		    ch[i-1+(k+jc*l1)*ido]=c1[i-1+(k+j*l1)*ido]+c1[i+(k+jc*l1)*ido];
		    ch[i+(k+j*l1)*ido]=c1[i+(k+j*l1)*ido]+c1[i-1+(k+jc*l1)*ido];
		    ch[i+(k+jc*l1)*ido]=c1[i+(k+j*l1)*ido]-c1[i-1+(k+jc*l1)*ido];
		}
	    }
	}
    }
    for(ik=0; ik<idl1; ik++)
	c2[ik]=ch2[ik];
    for(j=1; j<ip; j++)
	for(k=0; k<l1; k++)
	    c1[(k+j*l1)*ido]=ch[(k+j*l1)*ido];
    if(nbd<=l1)
    {
	is=-ido;
	for(j=1; j<ip; j++)
	{
	    is+=ido;
	    idij=is-1;
	    for(i=2; i<ido; i+=2)
	    {
		idij+=2;
		for(k=0; k<l1; k++)
		{
		    c1[i-1+(k+j*l1)*ido]=wa[idij-1]*ch[i-1+(k+j*l1)*ido]-wa[idij]*ch[i+(k+j*l1)*ido];
		    c1[i+(k+j*l1)*ido]=wa[idij-1]*ch[i+(k+j*l1)*ido]+wa[idij]*ch[i-1+(k+j*l1)*ido];
		}
	    }
	}
    }
    else
    {
	is=-ido;
	for(j=1; j<ip; j++)
	{
	    is+=ido;
	    for(k=0; k<l1; k++)
	    {
		idij=is;
		for(i=2; i<ido; i+=2)
		{
		    idij+=2;
		    c1[i-1+(k+j*l1)*ido]=wa[idij-1]*ch[i-1+(k+j*l1)*ido]-wa[idij]*ch[i+(k+j*l1)*ido];
		    c1[i+(k+j*l1)*ido]=wa[idij-1]*ch[i+(k+j*l1)*ido]+wa[idij]*ch[i-1+(k+j*l1)*ido];
		}
	    }
	}
    }
}				/*radbg*/

/*----------------------------------------------------------------------
   cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
  ----------------------------------------------------------------------*/

static void 
cfftf1(int n, double c[], double ch[], const double wa[], const int ifac[], int isign)
{
    int     idot, i;
    int     k1, l1, l2;
    int     na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;

    nf=ifac[1];
    na=0;
    l1=1;
    iw=0;
    for(k1=2; k1<=nf+1; k1++)
    {
	ip=ifac[k1];
	l2=ip*l1;
	ido=n / l2;
	idot=ido+ido;
	idl1=idot*l1;
	if(ip==4)
	{
	    ix2=iw+idot;
	    ix3=ix2+idot;
	    if(na==0)
		passf4(idot, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], isign);
	    else
		passf4(idot, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], isign);
	    na=1-na;
	}
	else if(ip==2)
	{
	    if(na==0)
		passf2(idot, l1, c, ch, &wa[iw], isign);
	    else
		passf2(idot, l1, ch, c, &wa[iw], isign);
	    na=1-na;
	}
	else if(ip==3)
	{
	    ix2=iw+idot;
	    if(na==0)
		passf3(idot, l1, c, ch, &wa[iw], &wa[ix2], isign);
	    else
		passf3(idot, l1, ch, c, &wa[iw], &wa[ix2], isign);
	    na=1-na;
	}
	else if(ip==5)
	{
	    ix2=iw+idot;
	    ix3=ix2+idot;
	    ix4=ix3+idot;
	    if(na==0)
		passf5(idot, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
	    else
		passf5(idot, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
	    na=1-na;
	}
	else
	{
	    if(na==0)
		passf(&nac, idot, ip, l1, idl1, c, c, c, ch, ch, &wa[iw], isign);
	    else
		passf(&nac, idot, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw], isign);
	    if(nac !=0)
		na=1-na;
	}
	l1=l2;
	iw+=(ip-1)*idot;
    }
    if(na==0)
	return;
    for(i=0; i<2*n; i++)
	c[i]=ch[i];
}				/*cfftf1*/


void 
cfftf(int n, double c[], double wsave[])
{
    int     iw1, iw2;

    if(n==1)
	return;
    iw1=2*n;
    iw2=iw1+2*n;
    cfftf1(n, c, wsave, wsave+iw1,(int*)(wsave+iw2),-1);
}				/*cfftf*/


void 
cfftb(int n, double c[], double wsave[])
{
    int     iw1, iw2;

    if(n==1)
	return;
    iw1=2*n;
    iw2=iw1+2*n;
    cfftf1(n, c, wsave, wsave+iw1,(int*)(wsave+iw2),+1);
}				/*cfftb*/


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

    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;
    i=1;
    l1=1;
    for(k1=1; k1<=nf; k1++)
    {
	ip=ifac[k1+1];
	ld=0;
	l2=l1*ip;
	ido=n / l2;
	idot=ido+ido+2;
	ipm=ip-1;
	for(j=1; j<=ipm; j++)
	{
	    i1=i;
	    wa[i-1]=1;
	    wa[i]=0;
	    ld+=l1;
	    fi=0;
	    argld=ld*argh;
	    for(ii=4; ii<=idot; ii+=2)
	    {
		i+=2;
		fi+=1;
		arg=fi*argld;
		wa[i-1]=cos(arg);
		wa[i]=sin(arg);
	    }
	    if(ip>5)
	    {
		wa[i1-1]=wa[i-1];
		wa[i1]=wa[i];
	    }
	}
	l1=l2;
    }
}				/*cffti1*/


void 
cffti(int n, double wsave[])
{
    int     iw1, iw2;

    if(n==1)
	return;
    iw1=2*n;
    iw2=iw1+2*n;
    cffti1(n, wsave+iw1,(int*)(wsave+iw2));
}				/*cffti*/

/*----------------------------------------------------------------------
   rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Real FFTs.
  ----------------------------------------------------------------------*/

static void 
rfftf1(int n, double c[], double ch[], const double wa[], const int ifac[])
{
    int     i;
    int     k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;

    nf=ifac[1];
    na=1;
    l2=n;
    iw=n-1;
    for(k1=1; k1<=nf;++k1)
    {
	kh=nf-k1;
	ip=ifac[kh+2];
	l1=l2 / ip;
	ido=n / l2;
	idl1=ido*l1;
	iw-=(ip-1)*ido;
	na=1-na;
	if(ip==4)
	{
	    ix2=iw+ido;
	    ix3=ix2+ido;
	    if(na==0)
		radf4(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
	    else
		radf4(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
	}
	else if(ip==2)
	{
	    if(na==0)
		radf2(ido, l1, c, ch, &wa[iw]);
	    else
		radf2(ido, l1, ch, c, &wa[iw]);
	}
	else if(ip==3)
	{
	    ix2=iw+ido;
	    if(na==0)
		radf3(ido, l1, c, ch, &wa[iw], &wa[ix2]);
	    else
		radf3(ido, l1, ch, c, &wa[iw], &wa[ix2]);
	}
	else if(ip==5)
	{
	    ix2=iw+ido;
	    ix3=ix2+ido;
	    ix4=ix3+ido;
	    if(na==0)
		radf5(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
	    else
		radf5(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
	}
	else
	{
	    if(ido==1)
		na=1-na;
	    if(na==0)
	    {
		radfg(ido, ip, l1, idl1, c, c, c, ch, ch, &wa[iw]);
		na=1;
	    }
	    else
	    {
		radfg(ido, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw]);
		na=0;
	    }
	}
	l2=l1;
    }
    if(na==1)
	return;
    for(i=0; i<n; i++)
	c[i]=ch[i];
}				/*rfftf1*/


void 
rfftb1(int n, double c[], double ch[], const double wa[], const int ifac[])
{
    int     i;
    int     k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;