www.pudn.com > 医学算法.rar > gen.c
#include
#include
/*
* gen - generate simulated data for CT scanners
*
* Carl Crawford
* David Nahamoo
*
* Purdue University
* West Lafayette, Indiana 47907
*
* Feburary 26, 1980
*/
/*
syntax: gen [options]
VARIABLES: DEFAULT
r1: distance from origin to translate axis 2.0
r2: distance from origin to detector axis 2.0
r: radius of scan circle 1.0
gmma: half beam width of fan 6.0
chi: full beam width of fan 12.0
n: number of detectors (maximum 100) 10
if: input file name for ellipse data gen.d
of: output file name for003.dat
m: number of ellipses (maximum 10) all
rf: radon space file name (invokes -r) radon
nfl: number of flashes in translate 100
k: number of projections 1
ver: version of generation routine 2
offset: Detector offset (added 7/12/86 by Malcolm) 0.0
FLAGS: DEFAULT
-a: equally angular detectors off
-e: equally spaced detectors on
-p: don't print program variables off
-r: save radon space data off
-1: generate first generation data on
-2: generate second generation data off
-3: generate third generation data off
-5: generate fifth generation data off
-s: sort flash data into projection data off
-d: -s with second and fifth generation deletion off
-w: weight function in fifth generation off
*/
#define ME 10 /* max number of ellipses */
#define MD 1500 /* max number of detectors */
float r1 = 2.0; /* distance from origin to translate axis */
float r2 = 2.0; /* distance from origin to detector axis */
float r = 1.0; /* radius of scan circle */
float gmma = 6.0; /* half beam width of fan (in degrees) */
float offset = 0.0; /* Detector Offset (in detector widths) */
int n = 10; /* number of detectors */
char *ifn = "gen.d"; /* ellipse file */
char *ofn = "for003.dat"; /* output file */
int m = 1; /* number of ellipses */
int ms; /* 1=m set */
char *rfn = "radon"; /* radon file */
int space; /* 1=equal angular data */
int print = 1; /* 1=print program variables */
int radon; /* save radon information */
int gen = 1; /* generation of data */
float x[ME]; /* x intercept of ellipse */
float y[ME]; /* y intercept of ellipse */
float theta[ME]; /* rotation of ellipse */
float c[ME]; /* attenuation value of ellipse */
float a[ME]; /* semi-major axis */
float b[ME]; /* semi-minor axis */
float abc2[ME]; /* a * b * c * 2 */
float cost[ME]; /* cos(theta[i]) */
float sint[ME]; /* sin(theta[i]) */
int sort; /* transform flash data to projection data */
float pi = 3.141592654; /* pie in the face (banana cream) ! */
FILE *ifd; /* input file descriptor */
FILE *ofd; /* output file descriptor */
FILE *rfd; /* radon file descriptor */
float beta[MD]; /* angle of each detector */
float tanb[MD]; /* tan[beta[i]] */
float cosb[MD]; /* cos[beta[i]] */
float sinb[MD]; /* sin[beta[i]] */
float rdt[MD]; /* radon coordinates */
float rda[MD]; /* radon coordinates */
float flash[MD]; /* integral flash data */
float *sbp,*sb; /* pointers to flash buffer */
int endt; /* indicate end of projection */
int nfl = 100; /* number of flashes in translate */
float tau; /* position of source on translate axis */
float rho; /* angle of translate axis */
float drho; /* angle between projections */
float dtau; /* step between projection samples */
int k = 1; /* number of projections */
int nc; /* ray number */
int kc; /* projection number */
int nfla; /* actual nfl for second gen. */
float taus; /* start of translation */
int delete; /* delete second generation data at ends */
float rhos; /* angle to last detector */
float tausg; /* taus / sin(gmma) */
float omegat; /* frequency of translate axis in gen 5 */
int ver = 2; /* routine version */
int weight; /* 1=weight ver 2 fifth data */
float aa; /* scale factor for -w option */
char *table[] = {
"r1",
"r2",
"r",
"gamma",
"n",
"if",
"of",
"m",
"rf",
"nfl",
"k",
"chi",
"ver",
"offset",
0
};
char *gent[] = {
"first",
"second",
"third",
"fourth",
"fifth",
0
};
main(argc,argv)
int argc;
char **argv;
{
register i;
char cc;
float dgmma,alpha,dt,t;
argc--;
argv++;
while(argc){
if(**argv == '-')while(cc = *++*argv)switch(cc){
case 'a': /* equal angle data */
space = 1;
break;
case 'e': /* equal space data */
space = 0;
break;
case 'p': /* don't print program variables */
print = 0;
break;
case 'r': /* save radon space info */
radon = 1;
break;
case 's': /* sort flash data into proj */
sort = 1;
break;
case '2': /* second generation data */
gen = 2;
break;
case '1': /* first generation data */
gen = 1;
break;
case '3': /* third generation data */
gen = 3;
break;
case '5': /* fifth generation data */
gen = 5;
break;
case 'd': /* sort with deletion */
delete = 1;
sort = 1;
break;
case 'w': /* weight ver 2 fifth */
ver = 2;
weight = 1;
gen = 5;
break;
default:
fprintf(stderr,"bad flag: -%c\n",cc);
exit(1);
}else switch(comm(*argv)){
case 1: /* r1 */
r1 = atof(*argv + 3);
break;
case 2: /* r2 */
r2 = atof(*argv + 3);
break;
case 3: /* r */
r = atof(*argv + 2);
break;
case 4: /* gmma */
gmma = atof(*argv + 6);
break;
case 5: /* n */
n = atoi(*argv + 2);
break;
case 6: /* if */
ifn = *argv + 3;
break;
case 7: /* of */
ofn = *argv + 3;
break;
case 8: /* m */
ms = 1;
m = atoi(*argv + 2);
break;
case 9: /* rf */
rfn = *argv + 3;
radon = 1;
break;
case 10: /* nfl */
nfl = atoi(*argv + 4);
break;
case 11: /* k */
k = atoi(*argv + 2);
break;
case 12: /* chi */
gmma = 0.5 * atof(*argv + 4);
break;
case 13: /* ver */
ver = atoi(*argv + 4);
break;
case 14: /* offset */
offset = atof(*argv+7);
break;
}
argc--;
argv++;
}
if(n > MD || n < 1){
fprintf(stderr,"bad detector count\n");
exit(1);
}
if((ifd = fopen(ifn,"r")) == NULL){
fprintf(stderr,"can't open: %s\n",ifn);
exit(1);
}
for(i=0;i ME){
fprintf(stderr,"too many ellipses\n");
exit(1);
}
}
fclose(ifd);
if((ofd = fopen(ofn,"w")) == NULL){
fprintf(stderr,"can't create: %s\n",ofn);
exit(1);
}
switch(gen){ /* init things specific to generation of data */
case 1: /* first generation */
drho = pi / k;
nfl = n;
dtau = (r + r) / nfl;
tau = -r - dtau * 0.5 + offset*dtau;
n = 1;
break;
case 2: /* second generation */
dtau = (r + r) / nfl;
space = 1;
k = 90.0 / gmma;
drho = 2.0 * gmma * pi / 180.0;
rhos = rho = gmma * pi / 180.0 * (1.0 - 1.0 / n);
tau = taus = -(r + r1*sin(rho)) / cos(rho) + dtau * 0.5 +
offset*dtau;
nfla = -taus / r * nfl;
tau -= dtau;
break;
case 3: /* third generation */
drho = 2.0 * pi / k;
tau = 0.0;
rho = -drho;
r1 = r / tan(gmma * pi / 180.0);
break;
case 5: /* fifth generation */
switch(ver){
case 1: /* small angle */
dtau = (r + r) / nfl;
space = 1;
k = 180.0 / gmma;
rhos = rho = gmma * pi / 180.0 * (1.0 - 1.0 / n);
tau = taus = -(r + r1*sin(rho)) / cos(rho) + dtau * 0.5;
nfla = -taus / r * nfl;
tau -= dtau;
tausg = taus / sin(rhos);
omegat = 2.0 * rhos / (nfla - 1.0);
rho -= omegat;
break;
case 2: /* w/o filter */
gmma = 225. / k;
taus = 1.5 * r / cos(gmma * pi / 180.0);
r1 = 0.5 * r / sin(gmma*pi/180.0);
aa = 5.0 * taus / (3.0 * gmma * pi / 180.0) - r1;
dtau = 2.0 * taus / (nfl - 1);
omegat = 2.0 * gmma / 5.0 * pi / 180.0;
drho = 3.0 * omegat / nfl;
rho = -drho * 0.5;
taus = -taus;
tau = taus - dtau;
break;
default:
fprintf(stderr,"bad version\n");
exit(1);
}
break;
}
if(radon){
if((rfd = fopen(rfn,"w")) == NULL){
fprintf(stderr,"can't create: %s\n",rfn);
exit(1);
}
}
if(sort && !(gen == 2 || gen == 5)){
fprintf(stderr,"-s flag not valid in this generation\n");
exit(1);
}
if(sort){
if((sb = (float *)calloc(nfla * n,sizeof(float))) == NULL){
fprintf(stderr,"can't allocate space for sort buffer\n");
exit(1);
}
}
if(print){
printf("\tGEN DATA\n\n");
printf("r1 = %f\n",r1);
printf("r2 = %f\n",r2);
printf("r = %f\n",r);
printf("gmma = %f\n",gmma);
printf("offset = %f\n",offset);
printf("n = %d\n",n);
printf("m = %d\n",m);
printf("input file: %s\n",ifn);
printf("output file: %s\n",ofn);
if(radon)printf("radon file: %s\n",rfn);
printf("equally %s data\n",space? "angular" : "spaced");
printf("%s format data\n",sort? "projection" : "flash");
printf("%s generation data\n",gent[gen -1]);
if(gen == 5)printf("version %d\n",ver);
if(weight)printf("post weighted data\n");
printf("k = %d\n",k*((gen==2 || (gen==5 && ver == 1))?n : 1));
if(gen <= 2 || gen==5){
printf("nfl = %d\n",nfl);
if(gen==2 || (gen==5 && ver == 1))printf("nfla = %d\n",nfla);
}
printf("\n\tELLIPSE DATA\n");
printf("\n x y theta c a b\n");
for(i=0;i