tomo3d.c

射线追踪程序

  if (!f) {
     perror(NULL);
     return(errno);
  }
  fwrite(&N, sizeof(int), 1, f);
  fwrite(&DS, sizeof(float), 1, f);
  fflush(f);
  for (i=0; i<N; i++) {
      fwrite(&(path[i]),sizeof(struct path_t),1,f);
      fwrite(path[i].x,sizeof(float),path[i].nseg,f);
      fwrite(path[i].y,sizeof(float),path[i].nseg,f);
      fwrite(path[i].z,sizeof(float),path[i].nseg,f);
      fwrite(path[i].T,sizeof(float),path[i].nseg,f);
      fflush(f);
  }
  fclose(f);
  if (buf) free(buf);
  return(0);
}

int blocknum(float x, float y, float z, struct block_t *blk, int nbb, int test)
{
        int i,j,k;
        struct block_t *b;
        if (test>-1) {
                k=test;
                b=&(blk[k]);
                if ((z>=b->z)&&(z<=b->z+b->dz)&&(x>=b->x)&&(x<=b->x+b->dx)
                                        &&(y>=b->y)&&(y<=b->y+b->dy)) {
                                return(k);
                }
        } else test=0;

        if (test<nbb/2) {
                j=test;
        } else {
                j=(nbb-test);
        }
        for (i=0; i<j; i++) {
                k=test+i;
                b=&(blk[k]);
                if ((z>=b->z)&&(z<=b->z+b->dz)&&(x>=b->x)&&(x<=b->x+b->dx)
                                        &&(y>=b->y)&&(y<=b->y+b->dy)) {
                                return(k);
                }
                k=test-i;
                b=&(blk[k]);
                if ((z>=b->z)&&(z<=b->z+b->dz)&&(x>=b->x)&&(x<=b->x+b->dx)
                                        &&(y>=b->y)&&(y<=b->y+b->dy)) {
                                return(k);
                }

                k=test+guess_blocknum_shift+i;
                if ((k>0)&&(k<nbb)) {
                        b=&(blk[k]);
                        if ((z>=b->z)&&(z<=b->z+b->dz)&&(x>=b->x)&&(x<=b->x+b->dx)
                                        &&(y>=b->y)&&(y<=b->y+b->dy)) {
                                return(k);
                        }
                }
                k=test+guess_blocknum_shift-i;
                if ((k>0)&&(k<nbb)) {
                        b=&(blk[k]);
                        if ((z>=b->z)&&(z<=b->z+b->dz)&&(x>=b->x)&&(x<=b->x+b->dx)
                                        &&(y>=b->y)&&(y<=b->y+b->dy)) {
                                return(k);
                        }
                }
        }


        if (test<nbb/2) {
                i=test+i;
                j=nbb;
        } else {
                j=test-i+1;
                i=0;
        }

        for (; i<j; i++) {
                b=&(blk[i]);
                if ((z>=b->z)&&(z<=b->z+b->dz)&&(x>=b->x)&&(x<=b->x+b->dx)
                                &&(y>=b->y)&&(y<=b->y+b->dy)) {
                        return(i);
                }
        }
        return(-1);
}

struct mod3d_t *readveloc(char *s, float dd, float W)
{
        struct mod3d_t *m3d=NULL;
        FILE *fxyz;
        int i,j,k,ii=0;
        struct mod_t *mod,*tmod;
        float minx,maxx,miny,maxy,maxz;
        m3d=(struct mod3d_t*)malloc(sizeof(struct mod3d_t));
        m3d->d=dd;
        m3d->W=W;
        fprintf(stderr,"Reading %s...",s);
        if ((fxyz=fopen(s,"rt"))==NULL) {
                perror(s);
                exit(0);
        }
        tmod=mod=(struct mod_t*)calloc(1,sizeof(struct mod_t));
        m3d->inmod=tmod;
        fscanf(fxyz,"%f %f %f %f %*s %*s %*s %*s",
                        &(mod->X),&(mod->Y),&(mod->Z),&(mod->V));
        while (!feof(fxyz)) {
                tmod->next=(struct mod_t*)calloc(1,sizeof(struct mod_t));
                if (fscanf(fxyz,"%f %f %f %f %*s %*s %*s %*s",
                                        &(tmod->next->X),&(tmod->next->Y),
                                        &(tmod->next->Z),&(tmod->next->V))!=4) {
                        free(tmod->next);
                        tmod->next=NULL;
                        break;
                }
                tmod=tmod->next;
        }

        /* back to absolute velocity */
        i=0;
        minx=maxx=mod->X;
        miny=maxy=mod->Y;
        maxz=mod->Z;
        for (tmod=mod; tmod; tmod=tmod->next) {
                tmod->V=Pvelocity(tmod->Z)*(1.0+0.01*tmod->V);

                /* ICI */
                /*
                 * tmod->V=Pvelocity(tmod->Z)*(1.0+0.1*tmod->V);
                 * */

                if (tmod->X<minx) minx=tmod->X;
                if (tmod->X>maxx) maxx=tmod->X;

                if (tmod->Y<miny) miny=tmod->Y;
                if (tmod->Y>maxy) maxy=tmod->Y;

                if (tmod->Z>maxz) maxz=tmod->Z;
                i++;
        }
        fprintf(stderr,"%d\n",i);
        m3d->X0=minx;
        m3d->Y0=miny;
        m3d->Nx=1+(maxx-minx)/dd;
        m3d->Ny=1+(maxy-miny)/dd;
        m3d->Nz=1+maxz/dd;
        m3d->mod=(float***)malloc(m3d->Nz*sizeof(float **));
        m3d->Nmod=(int***)malloc(m3d->Nz*sizeof(int **));
        for (i=0; i<m3d->Nz; i++) {
                m3d->mod[i]=(float**)malloc(m3d->Nx*sizeof(float *));
                if (!m3d->mod[i]) {
                        fprintf(stderr,"malloc error (readveloc %d\n",i);
                        exit(1);
                }
                m3d->Nmod[i]=(int**)malloc(m3d->Nx*sizeof(int *));
                if (!m3d->Nmod[i]) {
                        fprintf(stderr,"malloc error (readveloc %d\n",i);
                        exit(1);
                }
                for (j=0; j<m3d->Nx; j++) {
                        m3d->mod[i][j]=(float*)malloc(m3d->Ny*sizeof(float));
                        if (!m3d->mod[i][j]) {
                                fprintf(stderr,"malloc error (readveloc %d:%d\n",i,j);
                                exit(1);
                        }
                }
                for (j=0; j<m3d->Nx; j++) {
                        m3d->Nmod[i][j]=(int*)malloc(m3d->Ny*sizeof(int));
                        if (!m3d->Nmod[i][j]) {
                                fprintf(stderr,"malloc error (readveloc %d:%d\n",i,j);
                                exit(1);
                        }
                }
        }
        fprintf(stderr,"GRID: Z:%d X:%d Y:%d -> %d\n",
                        m3d->Nz,m3d->Nx,m3d->Ny,
                        m3d->Nz*m3d->Nx*m3d->Ny);

        fprintf(stderr,"Resampling...");
        for (i=0; i<m3d->Nz; i++) {
                for (j=0; j<m3d->Nx; j++) {
                        for (k=0; k<m3d->Ny; k++) {
                                m3d->mod[i][j][k]=0.0;
                        }
                }
        }

        for (tmod=m3d->inmod; tmod; tmod=tmod->next) {
                i=(int)(tmod->Z/m3d->d);
                j=(int)((tmod->X-m3d->X0)/m3d->d);
                k=(int)((tmod->Y-m3d->Y0)/m3d->d);
                if ((i<0)||(j<0)||(k<0)||
                                (i>=m3d->Nz)||(j>=m3d->Nx)||(k>=m3d->Ny)) {
                        fprintf(stderr,"oups: i=%d/%d j=%d/%d k=%d/%d\n",
                                        i,m3d->Nz,m3d->Nx,j,k,m3d->Ny);
                }
                m3d->Nmod[i][j][k]++;
                m3d->mod[i][j][k]+=tmod->V;
        }
        for (i=0; i<m3d->Nz; i++) {
                for (j=0; j<m3d->Nx; j++) {
                        for (k=0; k<m3d->Ny; k++) {
                                if (m3d->Nmod[i][j][k]) {
                                        m3d->mod[i][j][k]/=(float)m3d->Nmod[i][j][k];
                                        ii++;
                                } else {
                                        m3d->mod[i][j][k]=-1000.0;
                                }
                        }
                }
        }
        fprintf(stderr,"%d/%d\n",ii,m3d->Nz*m3d->Nx*m3d->Ny);

        if (i<5) {
                m3d->mod=NULL;
        }
        return(m3d);
}

float getvelocity(float x, float y, float z, struct mod3d_t *m3d, int test)
{
        int i,j,k;
        if (!m3d) {
                return(Pvelocity(z));
        }
        i=(int)(z/m3d->d);
        if (i>=m3d->Nz) {
                return(Pvelocity(z));
        }
        j=(int)((x-m3d->X0)/m3d->d);
        k=(int)((y-m3d->Y0)/m3d->d);
        if ((i<0)||(j<0)||(k<0)||
                        (i>=m3d->Nz)||(j>=m3d->Nx)||(k>=m3d->Ny)) {
                /*fprintf(stderr,"oups: i=%d j=%d k=%d\n",i,j,k);*/
                return(Pvelocity(z));
        }

        if (!m3d->Nmod[i][j][k]) {
                m3d->mod[i][j][k]=smoothvel(x,y,z,m3d,test);
                m3d->Nmod[i][j][k]++;
        }

        return(m3d->mod[i][j][k]);
}

float smoothvel(float x,float y,float z,struct mod3d_t *m3d,int test)
{
        double dist,vel,tdist,tmp,W2;
        struct mod_t *mod=m3d->inmod;
        struct mod_t *tmod;
        float X=1.0;
        W2=-0.25*(m3d->W)*(m3d->W)/log(0.5);
        do {
                dist=tdist=vel=0.;
                for (tmod=mod; tmod; tmod=tmod->next) {
                        dist=0.0;

                        tmp=(tmod->Z-z);
                        dist+=tmp*tmp;
                        if (dist>X*W2) continue;

                        tmp=(tmod->Y-y);
                        dist+=tmp*tmp;
                        if (dist>X*W2) continue;

                        tmp=(tmod->X-x);
                        dist+=tmp*tmp;
                        if (dist>X*W2) continue;

                        dist=exp(-0.5*dist/W2);
                        vel+=dist*tmod->V;
                        tdist+=dist;
                }
                X*=2.0;
                if (X>10.) {
                        float V=Pvelocity(z);
                        vel/=tdist;
                        if (fabs(vel/V)<0.15) {
                                return(vel);
                        } else {
                                return(V);
                        }
                }
        } while (tdist<EPS);
        return(vel/tdist);
}