www.pudn.com > 医学算法.rar > filt.c
#include
#include
/*
filt - filter tomographic projection data
Carl Crawford
Purdue University
West Lafayette, IN. 47907
March 17, 1980
*/
/*
Syntax: filt [options]
VARIABLES DEFAULT
k: number of projections 1
n: number of samples per projection 100
chi: full fan width 45.0
gamma: half fan width 22.5
m: log base 2 of NFFT 8
if: input file for003.dat
of: output file for002.dat
t: half duration of projection 1.0
r: radius of scan circle 1.0
nfl: number of flashes 100
FLAGS DEFAULT
-h: hamming window off
-k: kak's radial algorithm off
-x: lax's radial algorithm off
-e: equal spaced data on
-a: equal angular data off
-p: don't print program information off
-t: time domain filter on
-f: frequency domain filter off
-5: fifth generation data off
-s: save filter in 'filter' off
*/
int k = 1; /* number of projections */
int n = 100; /* rays per projection */
float chi = 45.0; /* full fan width */
int m = 8; /* log base 2 of NFFT */
int nfft; /* number of points in fft */
char *ifn = "for003.dat"; /* input file */
char *ofn = "for002.dat"; /* output file */
float t = 1.0; /* half duration of projections */
float r = 1.0; /* radius of scan circle */
int np; /* related to n */
int ham; /* 1 = use hamming window */
int nfl = 100; /* number of flashes */
int algo; /* 0=par 1=kak 2=lax 5=fifth */
int space; /* 1 = equal angular data */
int time; /* 1 = frequency domain filter */
int print = 1; /* 1 = print program info */
FILE *ifd; /* input file descriptor */
FILE *ofd; /* output file descriptor */
int odd; /* 1 = odd k */
float pi = 3.141592654; /* in indiana, pi = 3.0 */
float *pr1; /* pointer to projection one */
float *pr2; /* pointer to projection two */
float *wfilt; /* pointer to filter */
float *wgh; /* pointer to weight function */
int sfilt; /* 1 =save filter in 'filter' */
FILE *ffd; /* filter file descriptor */
char *table[] = {
"k",
"n",
"chi",
"gamma",
"m",
"if",
"of",
"t",
"r",
"nfl",
0
};
main(argc,argv)
int argc;
char **argv;
{
char cc;
register i,j;
float da,dt,dw,x,theta,rsc,tmp;
while(++argv,--argc){
if(**argv == '-')while(cc = *++*argv)switch(cc){
case 'h': /* hamming window */
ham = 1;
break;
case 'k': /* kak's algorithm */
algo = 1;
break;
case 'x': /* lax's algorithm */
algo = 2;
break;
case 'e': /* equal spaced data */
space = 0;
break;
case 'a': /* equal angular data */
space = 1;
break;
case 'p': /* don't print program info */
print = 0;
break;
case 't': /* time domain filter */
time = 0;
break;
case 'f': /* frequency domain filter */
time = 1;
break;
case '5': /* fifth generation data */
algo = 5;
break;
case 's': /* save filter */
sfilt = 1;
break;
default:
fprintf(stderr,"bad flag -%c\n",cc);
exit(1);
}else switch(comm(*argv)){
case 1: /* k */
k = atoi(*argv + 2);
break;
case 2: /* n */
n = atoi(*argv + 2);
break;
case 3: /* chi */
chi = atof(*argv + 4);
break;
case 4: /* gamma */
chi = 2.0 * atof(*argv + 6);
break;
case 5: /* m */
m = atoi(*argv + 2);
break;
case 6: /* if */
ifn = *argv + 3;
break;
case 7: /* of */
ofn = *argv + 3;
break;
case 8: /* t */
t = atof(*argv + 2);
break;
case 9: /* r */
r = atof(*argv + 2);
break;
case 10: /* nfl */
nfl = atoi(*argv + 4);
algo = 5;
break;
}
}
if( (ifd = fopen(ifn,"r") ) == NULL){
fprintf(stderr,"can't open: %s\n",ifn);
exit(1);
}
if( (ofd = fopen(ofn,"w") ) == NULL){
fprintf(stderr,"can't create: %s\n",ofn);
exit(1);
}
if(k & 1)odd = 1;
np = n;
if(algo == 5){
chi = 450.0 / k;
tmp = chi * pi / 180.0;
n = (4.0* asin(sin(tmp)/sqrt(5.0-cos(tmp)*4.0))/tmp+1.0)*n+0.5;
if((n - np) & 1)n++;
chi = (chi * n) / np;
if(!space){
n = tan(chi * 0.5 * pi / 180.0) / tan(tmp * 0.5) * np;
if((n - np) & 1) n++;
}
}
if((nfft = pow(2.0,(double)m) + 0.5) < (2*n - 1)){
fprintf(stderr,"WARNING: NFFT < 2N-1\n\n");
}
if((pr1=(float*)calloc(3*nfft +n*(algo > 0 ? 1 : 0),sizeof(float))) == NULL){
fprintf(stderr,"can't allocate buffer space\n");
exit(1);
}
pr2 = pr1 + nfft;
wfilt = pr2 + nfft;
wgh = wfilt + nfft;
if(print){
printf("FILT DATA\n\n");
printf("k = %d (%s)\n",k,(odd)? "odd" : "even");
printf("n = %d\n",n);
printf("chi = %f\n",chi);
printf("t = %f\n",t);
printf("nfft = %d\n",nfft);
printf("input file: %s\n",ifn);
printf("output file: %s\n",ofn);
printf("equal %s data\n",(space)? "angular" : "spaced");
printf("%s domain filter\n",(time)? "frequency" : "time");
if(sfilt)printf("filter saved in: %s\n","filter");
switch(algo){
case 0:
printf("parallel algorithm\n");
break;
case 1:
printf("kak's algorithm\n");
break;
case 2:
printf("lax's algorithm\n");
break;
case 5:
printf("fifth generation algorithm\n");
printf("np = %d\n",np);
printf("nfl = %d\n",nfl);
break;
}
printf("\n");
}
chi *= pi / 180.0;
rsc = 1.0 / (sin(0.5*chi) / cos(0.5*chi));
dt = 2.0 * t / n;
if(!space){ /* set up weight function */
switch(algo){
case 0:
case 5:
break;
case 1:
case 2:
x = -t + 0.5 * dt;
for(i=0;i>1);i += 2){
if(space){
pr1[i] = -da/(pi*sin(i*da)*sin(i*da)*k);
}else{
pr1[i] = -1.0 / (i*i*pi*dt*k);
}
pr1[nfft - i] = pr1[i];
}
fft(pr1,pr2,m,0);
for(i=0;i= 2){
fprintf(stderr,"no frequecy domain filter available\n");
exit(1);
}
dw = pi * n / (k * nfft * 2.0);
wfilt[0] = 0.0;
for(i=1;i<=(nfft>>1);i++){
wfilt[nfft - i] = wfilt[i] = i * dw;
}
}
if(ham){
for(i=1;i<=(nfft>>1);i++){
wfilt[nfft - i] = wfilt[i] *= cos(pi*i/nfft);
}
}
if(sfilt){
if((ffd = fopen("filter","w")) == NULL){
fprintf(stderr,"can't create: %s\n","filter");
exit(1);
}
fwrite(wfilt,sizeof(float),nfft,ffd);
fclose(ffd);
exit(0);
}
if(algo == 5){
k *= nfl;
if(k & 1)odd = 1;
else odd = 0;
}
for(i=0;i<((k+odd)>>1);i++){
input(pr1);
if(i || !odd)input(pr2);
if(algo >= 1 && algo <=4){
for(j=0;j