www.pudn.com > 医学算法.rar > back.c
#include
#include
/*
back - back projector
Carl Crawford
Purdue University
W. Lafayette, IN. 47907
March 6, 1978
*/
/*
syntax: back [options]
VARIABLES DEFAULT
v: number of projections 50
k: same as v: 50
n: number of rays per projection 100
d: input file for002.dat
if: same as d: for002.dat
l: number of lines in -l option 1
lines: same as l: 1
skip: number of projections to ignore 0
x: line in space for -l 0.0
y: same as x: 0.0
t: initial object rotation 0.0
theta: same as t: 0.0
c: projection center offset 0.0
cen: same as c: 0.0
gx: x start for grid -1.0
gy: x start for grid 1.0
gs: width of grid 2.0
chi: full fan width 45.0
gamma: half fan width 22.5
q: quadrant of reconstuction entire
of: output file pic.r
r: radius of reconstruction circle 1.0
nfl: number of flashes for fifth generation 100
lf: output line in this file line
FLAGS DEFAULT
-x: lax's algorithm off
-k: kak's algorithm off
-o: assume circularily symmetrical object off
-l: reconstruct only a line off
-c: shift grid by half step off
-5: fifth generation data off
-p: don't print program status off
-a: equal angular data off
-e: equal space data on
-f: filtered fifth (invokes -5) off
*/
int algo; /* 0=parallel 1=kak 2=lax 5=fifth */
int one; /* 1=use single projection */
int line; /* 1=reconstruct only line of picture */
int views = 50; /* default number of views */
int rays = 100; /* default number of rays */
int rays1; /* rays - 1 */
char *data = "for002.dat"; /* default input file name */
char *ofn = "pic.r"; /* default output file name */
int npic = 64; /* size of grid, npic * npic */
int npic2 = 128; /* size of grid for line picture */
int space; /* 1 = equal angular data */
int print = 1; /* 1 = print program info */
float pic[4096],*p; /* picture and pointer to it */
float proj[512]; /* projection data */
float bb[512]; /* break points for equal angle data */
float theta,dtheta; /* angle value for each view */
float offset; /* projection offset */
float grid,grido2; /* size of grid block */
float pi = 3.1415926535;
float convert; /* conversion between degress and radians */
float cost, /* cos(theta) */
sint, /* sin(theta) */
cosd,cosd2, /* cost * grid */
sind,sind2; /* sint * grid */
int l; /* integer location of projection value */
int ll; /* related to l */
int fd; /* input data fd and output fd */
int ii1,ii2; /* set if in line routine */
float yy1s = 0.0; /* default value of line reconstruction */
float yy1; /* related to yy1s */
float gridx = -1.0; /* x starting grid position */
float gridy = 1.0; /* y starting grid position */
float grids = 2.0; /* size of the grid */
int lines = 1; /* number of line sections in the -l option */
int skip = 0; /* skip skip projections between reading */
int psf; /* 1=make center grid point 0,0 */
int quad; /* quadrant reconstruction */
char *lfn = "line"; /* line file name */
float radius = 1.0; /* radius of reconstruction circle */
/* variables used only in the kak and lax algorithm */
float rdc, /* distance between the center and the screen */
r1, /* distance from detectors to source */
rsc, /* distance fron source to the center */
half, /* set so that delta t in the proj. is one */
chi = 45.0, /* fan width */
t1, /* like t */
t1s, /* related to t1 */
t1r1, /* t * r1 */
t1r1s, /* modified t1r1 */
s1, /* s */
s1s, /* modified s1 */
cosdr1, /* cosd * r1 */
sindr1, /* sind * r1 */
cosdrsc, /* cosd / rsc */
sindrsc, /* sind / rsc */
s1rsc, /* s1 / rsc */
s1rscs, /* related to s1rsc */
tpr; /* location in the projection */
/* variables used only in the parallel algorithm */
float r, /* location in the projection */
rs, /* modified r */
t, /* delta in the projections */
gridt; /* grid / t */
/* variables used only for fifth generation */
int nfl = 100; /* number of flashes */
float drho; /* step between flashes */
float taus; /* start of translation */
float tau; /* translate position */
float dtau; /* step between flashes */
int filt; /* 1=filtered fifth */
int np; /* related to n */
char *table[] = {
"v",
"k",
"r",
"n",
"d",
"if",
"l",
"lines",
"skip",
"x",
"y",
"t",
"theta",
"c",
"cen",
"gx",
"gy",
"gs",
"chi",
"gamma",
"q",
"of",
"nfl",
"lf",
0
};
char *algon[] = {
"parallel",
"kak's radial",
"lax's radial",
"",
"",
"fifth generation"
};
main(argc,argv)
int argc;
char **argv;
{
char cc;
register int i,j;
while(++argv,--argc){
if(**argv == '-')while(cc = *++*argv)switch(cc){
case 'x': /* use lax's algorithm */
algo = 2;
break;
case 'k': /* use kak's algorithm */
algo = 1;
break;
case 'o': /* use single projection */
one = 1;
break;
case 'l': /* reconstruct only line of picture */
line = 1;
break;
case 'c': /* make center grid 0,0 */
psf = 1;
break;
case '5': /* fifth generation data */
algo = 5;
break;
case 'p': /* don't print program info */
print = 0;
break;
case 'a': /* equal angular data */
space = 1;
break;
case 'e': /* equal space data */
space = 0;
break;
case 'f': /* filtered fifth */
filt = 1;
algo = 5;
break;
default:
fprintf(stderr,"bad flag: -%c\n",cc);
exit(1);
}else switch(comm(*argv)){
case 1: /* v */
case 2: /* k */
views = atoi(*argv+2);
break;
case 3: /* r */
radius = atof(*argv + 2);
break;
case 4: /* n */
rays = atoi(*argv+2);
break;
case 5: /* d */
data = *argv + 2;
break;
case 6: /* if */
data = *argv + 3;
break;
case 7: /* l */
lines = atoi(*argv + 2);
line = 1; /* Set the -l option */
break;
case 8: /* lines */
lines = atoi(*argv + 5);
break;
case 9: /* skip */
skip = atoi(*argv + 5);
break;
case 10: /* x */
case 11: /* y */
yy1s = atof(*argv + 2);
break;
case 12: /* t */
theta = atof(*argv + 2);
break;
case 13: /* theta */
theta = atof(*argv + 6);
break;
case 14: /* c */
offset = atof(*argv + 2);
break;
case 15: /* cen */
offset = atof(*argv + 4);
break;
case 16: /* gx */
gridx = atof(*argv + 3);
break;
case 17: /* gy */
gridy = atof(*argv + 3);
break;
case 18: /* gs */
grids = atof(*argv + 3);
break;
case 19: /* chi */
chi = atof(*argv + 4);
break;
case 20: /* gamma */
chi = 2.0 * atof(*argv + 6);
break;
case 21: /* q */
quad = atoi(*argv + 2);
break;
case 22: /* of */
ofn = *argv + 3;
break;
case 23: /* nfl */
nfl = atoi(*argv + 4);
break;
case 24: /* lf */
lfn = *argv + 3;
line = 1;
break;
}
}
init();
if(line)ofn = lfn;
if(print){
printf("BACK DATA\n\n");
printf("n = %d\n",rays);
printf("k = %d\n",views);
printf("r = %f\n",radius);
printf("gx = %f\n",gridx);
printf("gy = %f\n",gridy);
printf("gs = %f\n",grids);
printf("theta = %f\n",theta/convert);
printf("cen = %f\n",offset);
printf("%s back projection algorithm\n",algon[algo]);
printf("input file: %s\n",data);
printf("output file: %s\n",ofn);
printf("equal %s data\n",(space)? "angular" : "spaced");
if(psf)printf("psf mode on\n");
if(line){
printf("line reconstruction: ");
printf("%d line(s)",lines);
printf(" starting at y = %f\n",yy1s);
}
if(one)printf("one projection mode on\n");
if(skip)printf("skip %d projection(s)\n",skip);
if(algo){
printf("chi = %f\n",chi/convert);
printf("rsc = %f\n",rsc);
printf("rdc = %f\n",rdc);
if(filt)printf("np = %d\n",np);
}
if(algo == 5){
printf("nfl = %d\n",nfl);
printf("%sfiltered version\n",(filt)? "" : "un");
}
printf("\n");
}
if(one){
readp();
for(i=0;i=0 && l =0 && l= rs){
while(r > bb[ll]){
if(ll == rays1)break;
ll++;
}
l = ll;
}else{
while(r < bb[ll]){
if(ll == 0)break;
ll--;
}
l = ll + 1;
}
if(l>1 && l=0 && l=0 && l=0 && l= rs){
while(r > bb[ll]){
if(ll == rays1)break;
ll++;
}
l = ll;
}else{
while(r < bb[ll]){
if(ll == 0)break;
ll--;
}
l = ll + 1;
}
if(l>1 && l