www.pudn.com > 医学算法.rar > path.c
#include
#include
#define NT 10 /* resolution of TOF waveforms */
#define SIZE 128 /* size of water path signal */
#define M 1024 /* length of tof buffer */
#define NM 128 /* size of projection buffer */
#define WL 256 /* input water length */
/*
path - generate multipath data
Carl Crawford
Purdue University
W. Lafayette, IN. 47907
Jan. 15, 1981
*/
float a = 4.0; /* semi - major axis (cm) */
float b = 4.0; /* semi - minor axis (cm) */
float width = 1.0; /* width of transducer (cm) */
float w2; /* width * 0.5 */
float r; /* same as width */
float h; /* used in area calculation */
float ar,ar2; /* areas */
int print; /* 1=don't print program status */
float pi = 3.14159265;/* pi in the face */
int n = 050; /* samples per projection */
int k = 1; /* number of projections */
float level = 2.; /* threshold level for tof */
char *ofn = "for003.dat"; /* output file name */
char *ifn = "water"; /* water path signal name */
float gam = 1.0; /* attenuation frequency dependence */
float gamp = 1.0; /* gamma prime */
float vw = 1500.; /* velocity of water (m/s) */
float c = 1550.; /* velocity of object (m/s) */
float alpha = .05; /* slope of attenuation */
FILE *ifd,*ofd; /* file descriptors */
FILE *sfd; /* specturm file */
FILE *fd1; /* for transform file */
int gen; /* 1= generate tof file */
float sigma = 3.0; /* standard deviation of power spectrum */
char *gfn = "tof"; /* tof file name */
int in[WL]; /* water path signal */
float *zr,*zi,*qr,*qi,*dh,*fg,*mag;
float dt = .01; /* delta time (usec) */
float df; /* delta freq. (mhz) */
float dis = 25.0; /* distance between transducers (cm) */
float dtau; /* bump between projection samples */
float tau; /* projection sample position */
float taus; /* first tau */
float proj[NM]; /* projection */
float theta; /* angular projection position */
float dtheta; /* angualr bump between projections */
float a2; /* a*a*cost*cost+b*b*sint*sint */
float aa; /* related to a2 */
float ab; /* 2 * a * b / a2 */
float ptt; /* raw projection value */
int time; /* time delay */
int sig; /* which tof attenuated signal */
int i1,i2; /* positioning for tof buffer */
float *f1,*f2; /* character pointers */
float time_scale; /* time relationship between tof buffer and time */
float *buf,*tbuf; /* buffer pointers */
int mult; /* 1=include multipath */
int fl1,fl2; /* indicate if transducer edges are within ellipse */
int tof; /* 1=generate tof data */
int att = 1; /* 1=generate attenuation data */
float m0,m1; /* moments */
float tmp; /* temp var */
int b1 = 3; /* band for division */
int b2 = 10; /* band for division */
int div; /* 1=spectral division */
float c0,c1; /* constants for spectral division */
int save; /* 1=save spectrum */
char *table[] = {
"a",
"b",
"width",
"n",
"level",
"gamma",
"vw",
"c",
"alpha",
"of",
"if",
"k",
"gf",
"dt",
"dis",
"gammap",
"sigma",
"b1",
"b2",
0
};
main(argc,argv)
int argc;
char **argv;
{
char cc;
register i;
while(argv++ , --argc){
if(argv[0][0] == '-')while(cc = *++*argv)switch(cc){
case 'p': /* print flag */
print = 1;
break;
case 'g': /* generate tof file */
gen = 1;
break;
case 'm': /* include multipath */
mult = 1;
break;
case 't': /* generate tof data */
tof = 1;
att = 0;
break;
case 'a': /* generate attenuation data */
att = 1;
tof = 0;
break;
case 'd': /* spectral division */
div = 1;
att = 1;
tof = 0;
break;
case 's': /* save spectrum */
save = 1;
break;
default:
fprintf(stderr,"bad flag: -%c\n",cc);
exit(1);
}
else
switch(i = comm(*argv)){
case 1:
a = atof(*argv + 2);
break;
case 2:
b = atof(*argv + 2);
break;
case 3:
width = atof(*argv + 6);
break;
case 4:
n = atoi(*argv + 2);
break;
case 5:
level = atof(*argv + 6);
break;
case 6:
gam = atof(*argv + 6);
break;
case 7:
vw = atof(*argv + 3);
break;
case 8:
c = atof(*argv + 2);
break;
case 9:
alpha = atof(*argv + 6);
break;
case 10:
ofn = *argv + 3;
break;
case 11:
ifn = *argv + 3;
break;
case 12:
k = atoi(*argv + 2);
break;
case 13:
gfn = *argv + 3;
break;
case 14:
dt = atof(*argv + 3);
break;
case 15:
dis = atof(*argv + 4);
break;
case 16:
gamp = atof(*argv + 7);
break;
case 17:
sigma = atof(*argv + 6);
break;
case 18:
b1 = atoi(*argv + 3);
break;
case 19:
b2 = atoi(*argv + 3);
break;
}
}
if(n <= 0 || n > NM)err("illegal n","");
ifn = gen? ifn : gfn;
ofn = gen? gfn : ofn;
if(!print){
printf("PATH DATA\n\n");
printf("a: %f (cm)\n",a);
printf("b: %f (cm)\n",b);
printf("width: %f (cm)\n",width);
printf("dis: %f (cm)\n",dis);
printf("n: %d\n",n);
printf("k: %d\n",k);
printf("level: %f\n",level);
printf("vw: %f (m/s)\n",vw);
printf("c: %f (m/s)\n",c);
printf("alpha: %f\n",alpha);
printf("sigma: %f\n",sigma);
printf("gamma: %f\n",gam);
printf("gammap: %f\n",gamp);
if(div)printf("b1: %d\n",b1);
if(div)printf("b2: %d\n",b2);
printf("dt: %f (usec)\n",dt);
printf("input file: %s\n",ifn);
printf("output file: %s\n",ofn);
if(gen)printf("generate tof file\n");
if(!gen)printf("generate %s data\n",tof? "tof":"attenuation");
if(att && !gen)printf("%s technique\n",div?"division" : "fs");
if(mult)printf("include multipath\n");
if(save)printf("spectrum saved in 'spectrum'\n");
}
df = 1.0 / (WL * dt);
w2 = width * 0.5;
r = w2;
if((ifd = fopen(ifn,"r")) == NULL)err("can't open: ",ifn);
if((ofd = fopen(ofn,"w")) == NULL)err("can't create: ",ofn);
if(save)if((sfd = fopen("spectrum","w")) == NULL)err("can't create: spectrum\n");
if(gen){
gtof();
}else{
gdat();
}
}
comm(s)
char *s;
{
register int i,j,r;
for(i=0;table[i];i++){
for(j=0;(r=table[i][j]) == s[j] && r;j++);
if(r == 0 && s[j] == '=' && s[j+1] )return(i+1);
}
fprintf(stderr,"bad option: %s\n",s);
exit(1);
}
err(s1,s2)
char *s1,*s2;
{
fprintf(stderr,"%s%s\n",s1,s2);
exit(1);
}
gtof()
{
register i,j;
if((zr = (float *)calloc(5 * WL,sizeof(float))) == NULL){
err("can't allocate space","");
}
zi = zr + WL;
qr = zi + WL;
qi = qr + WL;
dh = qi + WL;
if(fread(in,sizeof(int),WL,ifd) != WL)
err("unexpected EOF in: ",ifd);
for(i=0;i c){ /* water faster than object */
i1 = 0;
i2 = SIZE;
}else{ /* object faster than water */
i2 = M;
i1 = i2 - SIZE;
}
dtau = dis / n;
taus = -(dis + dtau) * 0.5;
dtheta = pi / k;
theta = - dtheta;
time_scale = 10000 / dt * (1.0 / c - 1.0 / vw);
for(ii=0;ii= level)break;
if(i == M){
fprintf(stderr,"threshold not met. n: %d\n",j);
exit(1);
}
proj[j] = i;
}else{
for(i=0;i= aa){
ptt = 0.0;
}else{
ptt = ab * sqrt(a2 - pos * pos);
}
if(tof){
time = ptt * time_scale;
sig = floor(ptt * NT / (2.0 * a));
if(sig >= NT)sig = NT - 1;
if(abs(time) > (M - SIZE))err("too much delay","");
f2 = tbuf + i1 + time;
for(i=i1,f1=buf+sig*SIZE;i