parser.y

蒙特卡罗模拟光子成像C语言版,代码简洁专业

  if (dv <= 0) 
    dv = 1.0;           /* default when no length specified */

  config->dl = dv;

  /* Back-fill any missing lengths (implicitly the largest one) */

  if (config->xstep <= 0.0) config->xstep = config->dl;
  if (config->ystep <= 0.0) config->ystep = config->dl;
  if (config->zstep <= 0.0) config->zstep = config->dl;

  /* Calculate min and max photon lengths from the min and max
     propagation times in vacuuo */

  if (checkconfig() != 0)         /* Sanity check the final configuration */
    return 2;

  return 0;                       /* Done */
}

void init_config(void)
{
  int i;
#if (YYDEBUG)
  extern int yydebug;

  yydebug = 1;
#endif

  /* Set default configuration */

  config->v0 = C_VACUUM;

  config->xstep = config->ystep = config->zstep = -1;

  config->segFile[0] = '\0';

  config->Ixmin = -1;  config->Ixmax = -1;
  config->Iymin = -1;  config->Iymax = -1;
  config->Izmin = -1;  config->Izmax = -1;

  config->Ntissue = 0;

  for (i = 0; i < MAXTISSUE; i++)
    {
      config->tmus[i] = MUS_AIR;
      config->tmua[i] = MUS_AIR;
      config->tg[i]   =   G_AIR;
      config->tn[i]   = IDX_AIR;
    }

  config->nDets = 0;

  for (i = 0; i < MAXDET; i++)
    {
      config->detLoc[i][0] = -1;
      config->detLoc[i][1] = -1;
      config->detLoc[i][2] = -1;

      config->detDir[i][0] = -1;
      config->detDir[i][1] = -1;
      config->detDir[i][2] = -1;

      config->detRad[i]    = -1;
      config->detNA[i]     = -1;
    }

  config->nSrcs = 0;

  for (i = 0; i < MAXDET; i++)
    {
      config->srcLoc[i][0] = -1;
      config->srcLoc[i][1] = -1;
      config->srcLoc[i][2] = -1;

      config->srcDir[i][0] = -1;
      config->srcDir[i][1] = -1;
      config->srcDir[i][2] = -1;

      config->srcRad[i]    = -1;
      config->srcNA[i]     = -1;
    }

  return;
}

void add_tissue(double mus, double mua, double tg, double tn)
{
  int n = ++config->Ntissue;

  config->tmus[n] = mus;
  config->tmua[n] = mua;
  config->tg[n]   = tg;
  config->tn[n]   = tn;

  return;
}

/* Add a new detector entry based on the supplied detector parameters */

void add_detector(double r[3], double d[3], double rad, double na)
{
  double nrm = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
  int      n = config->nDets;

#if (HAVE_ASSERT_H)
  assert((n >= 0) && (n < MAXDET));
#endif

  if ((r[0] == -998) && (r[1] == -998) && (r[2] == -998)) 
    {
      yyerror("Detector position not defined\n");
      exit(1);
    }
  else
    {
      config->detLoc[n][0] = r[0];
      config->detLoc[n][1] = r[1];
      config->detLoc[n][2] = r[2];
    }

  if (nrm > 0)
    {
      /* Turn direction into a unit vector */

      config->detDir[n][0] = d[0]/nrm;
      config->detDir[n][1] = d[1]/nrm;
      config->detDir[n][2] = d[2]/nrm;
    }
  else
    {
      config->detDir[n][0] = 0;
      config->detDir[n][1] = 0;
      config->detDir[n][2] = 0;
    }

  if (rad <= 0)
    {
      yyerror("Detector radius non-positive or not defined");
      exit(1);
    }
  else
    config->detRad[n] = rad;

  if (na < 0)
    {
      yyerror("Detector NA must be non-negative");
      exit(1);
    }
  else if (na > 1)
    {
      yyerror("Source NA too large");
      exit(1);
    }
  else
    config->detNA[n] = na;

  ++config->nDets;

  return;
}

void add_source(double r[3], double d[3], double rad, double na)
{
  double nrm = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
  int      n = config->nSrcs;

#if (HAVE_ASSERT_H)
  assert((n >= 0) && (n < MAXDET));
#endif

  if ((r[0] == -998) && (r[1] == -998) && (r[2] == -998)) 
    {
      yyerror("Source position not defined\n");
      exit(1);
    }
  else
    {
      config->srcLoc[n][0] = r[0];
      config->srcLoc[n][1] = r[1];
      config->srcLoc[n][2] = r[2];
    }

  if (nrm > 0)
    {
      /* Turn direction into a unit vector */

      config->srcDir[n][0] = d[0]/nrm;
      config->srcDir[n][1] = d[1]/nrm;
      config->srcDir[n][2] = d[2]/nrm;
    }
  else
    {
      yyerror("Source direction zeros or not defined\n");
      exit(1);
    }

  if (rad < 0)
    {
      yyerror("Source radius negative or not defined");
      exit(1);
    }
  else
    config->srcRad[n] = rad;

  if (na < 0)
    {
      yyerror("Source NA must be non-negative");
      exit(1);
    }
  else if (na > PI/2)
    {
      yyerror("Source NA too large");
      exit(1);
    }
  else
    config->srcNA[n] = na;

  ++config->nSrcs;

  return;
}

#define PARSE_ERROR(s) { printf("Input error: %s\n", (s)); return 1; }
#define PARSE_WARN(s)  { printf("Warning: %s\n", (s)); }

int checkconfig(void)
{
  int i;

  if ((config->frequency > 1.0e10) || (config->frequency < -1.0e10))
    printf("Warning: abnormally large frequency reqested (%f MHz)\n",
           config->frequency * 1e-6);

  if ((config->nTstep <= 0) || (config->stepT <= 0))
    PARSE_ERROR("Non-causal time steps not allowed");

  if ((config->xstep <= 0) || (config->ystep <= 0) || (config->zstep <= 0))
    PARSE_ERROR("xstep, ystep, and zstep must be strictly positive");

  if ((config->xstep != config->ystep) || (config->xstep != config->zstep))
    PARSE_ERROR("xstep == ystep == zstep required in this version");

  if ((config->Ixmin < 0) || (config->nIxstep > config->nxstep))
    PARSE_ERROR("Illegal image size (X dimension)");

  if ((config->Iymin < 0) || (config->nIystep > config->nystep))
    PARSE_ERROR("Illegal image size (Y dimension)");

  if ((config->Izmin < 0) || (config->nIzstep > config->nzstep))
    PARSE_ERROR("Illegal image size (Z dimension)");

  if (config->Ntissue >= MAXTISSUE)
    PARSE_ERROR("Too many tissue types");

  for (i = 1; i <= config->Ntissue; i++)
    {
      if (config->tmus[i] == 0)
        PARSE_WARN("Zero scattering coefficient");
      if (config->tmus[i] <  0)
        PARSE_ERROR("Non-physical scattering coefficient");
      if (config->tmua[i] <  0)
        PARSE_ERROR("Non-physical absorption coefficient");
      if ((config->tg[i] < -0.5) || (config->tg[i] > 1.0))
        PARSE_ERROR("Non-physical anisotropy coefficient");
      if (config->tn[i] <= 0)
        PARSE_ERROR("Non-physical index of refraction");

      if (config->tn[i] < 1.0)
        PARSE_WARN(
          "Index of refraction < 1, I hope you know what you're doing!\n");
    }

  if (config->nSrcs <= 0)
    PARSE_ERROR("No Source Specified");

  for (i = 0; i < config->nSrcs; i++)
    if ((config->srcLoc[i][0] < 0) ||
	(config->srcLoc[i][1] < 0) ||
	(config->srcLoc[i][2] < 0) ||
	(config->srcLoc[i][0] >= config->nxstep * config->xstep) ||
	(config->srcLoc[i][1] >= config->nystep * config->ystep) ||
	(config->srcLoc[i][2] >= config->nzstep * config->zstep))
      {
	PARSE_WARN("Source located outside of image volume");
	fprintf(stderr, "\tr = [ %f,%f,%f ]\n", config->srcLoc[i][0],
		config->srcLoc[i][1], config->srcLoc[i][2]);
      }

  for (i = 0; i < config->nDets; i++)
    if ((config->detLoc[i][0] < 0) ||
	(config->detLoc[i][1] < 0) ||
	(config->detLoc[i][2] < 0) ||
	(config->detLoc[i][0] >= config->nxstep * config->xstep) ||
	(config->detLoc[i][1] >= config->nystep * config->ystep) ||
	(config->detLoc[i][2] >= config->nzstep * config->zstep))
      {
	PARSE_WARN("Detector located outside of image volume");
	fprintf(stderr, "\tr = [ %f,%f,%f ]\n", config->detLoc[i][0],
		config->detLoc[i][1], config->detLoc[i][2]);
      }

  return 0;
}