commonmc.c

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

      *ky /= knorm;
      *kz /= knorm;

#if (PARANOID)
      assert(SQR(*kx) + SQR(*ky) + SQR(*kz) > 0.9999);
      assert(SQR(*kx) + SQR(*ky) + SQR(*kz) < 1.0001);
#endif
    }

  return 0;
}

/**************************************************************
  HANDLE PHASE ROTATION BETWEEN REAL AND IMAGINARY PART
**************************************************************/

/* Originally, this handled just the increment due to the current
 *  step.  That proved sufficiently non-robust against round-off
 *  errors to warrent recomputing the photon amplitude and phase
 *  every step of the simulation.  Since I was already computing
 *  an exp(), a sin(), and a cos(), this is not significantly
 *  slower than the old way, which is also nice.
 */

static void rotate_phase(TWOPT_T *u, TWOPT_T *v, 
                         const struct Config *syscfg, const float *lenTiss)
{
  double len, phi, mag;
  int i;

  len = 0.0;
  phi = 0.0;

  for (i = 1; i <= syscfg->Ntissue; i++)
    {
      len += lenTiss[i] * syscfg->tmua[i];		/* \mu l */
      phi += lenTiss[i] * (syscfg->tn[i] * syscfg->k0);	/* phase */
    }

  mag = exp(-len);

#if (PARANOID)
  assert(mag > 0);
#endif

  if (syscfg->frequency == 0.0)		/* CW case */
    {
      *u = mag;
      *v = 0;
    }
  else
    /* RF case */
    {
      *u = mag * cos(phi);
      *v = mag * sin(phi);
    }

  return;
}

/**************************************************************
 Renormalized coordinte structure to put fractional part back
  on to 0 <= f < 1.
**************************************************************/

static void renorm(MCCOORD *r)
{
  /* Move the fractional portions back to 0 <= f < 1. Adjust integer 
   * portions to match
   */

  if ((r->fx < 0.0) || (r->fx >= 1.0))
    {
      for ( ; r->fx  < 0.0; r->fx += 1.0) r->ix -= 1;
      for ( ; r->fx >= 1.0; r->fx -= 1.0) r->ix += 1;
    }

  if ((r->fy < 0.0) || (r->fy >= 1.0))
    {
      for ( ; r->fy  < 0.0; r->fy += 1.0) r->iy -= 1;
      for ( ; r->fy >= 1.0; r->fy -= 1.0) r->iy += 1;
    }

  if ((r->fz < 0.0) || (r->fz >= 1.0))
    {
      for ( ; r->fz  < 0.0; r->fz += 1.0) r->iz -= 1;
      for ( ; r->fz >= 1.0; r->fz -= 1.0) r->iz += 1;
    }

  return;
}

/**************************************************************
 *
 * CALCULATE THE NEW SCATTERING ANGLE USING HENYEY-GREENSTEIN
 *
 *************************************************************/

static void scatter(double gg, double *kx, double *ky, double *kz, int *seed)
{
#if (PARANOID)
  /* Make sure the anisotropy is physically meaningful */

  assert((gg >= -0.5) && (gg <= 1.0));

  /* Sanity-check the input direction cosines */

  assert((*kx >= -1) && (*kx <= 1));
  assert((*ky >= -1) && (*ky <= 1));
  assert((*kz >= -1) && (*kz <= 1));
  assert(SQR(*kx) + SQR(*ky) + SQR(*kz) > 0.9999);
  assert(SQR(*kx) + SQR(*ky) + SQR(*kz) < 1.0001);
#endif

  /* If gg is identically zero, the calculated distribution has a
   * divide-by-zero, but gg==0 is just the isotropic case, which is
   * trivial to handle */

  if (gg == 0.0)
    {
      double nrm, C1, C2, C3;

      /* Select in a cube, reject if it falls outside the sphere.
         Gives us a vector uniformly distributed through-out the
         volume of the sphere.  Rescale to unit length after.  May not
         be the most efficient algorithm imaginable, but it does
         guarantee uniform angular sampling. */

      do
	{
	  C1 = 2 * ran(0) - 1;
	  C2 = 2 * ran(0) - 1;
	  C3 = 2 * ran(0) - 1;
	}
      while (SQR(C1) + SQR(C2) + SQR(C3) > 1.0);

      /* (kx,ky,kz) must be a unit vector, rescale vector length */
      nrm = sqrt(SQR(C1) + SQR(C2) + SQR(C3));

#if (PARANOID)
      assert(nrm > 0.0);
#endif

      *kx = C1 / nrm;
      *ky = C2 / nrm;
      *kz = C3 / nrm;
    }
  else
    {			      /* gg != 0 */
      double kxo, kyo, kzo;   /* OLD DIRECTION COSINES */
      double phi, cphi, sphi, theta, stheta, ctheta;
      double rnm, nrm, foo;

      kxo = *kx;
      kyo = *ky;
      kzo = *kz;

      phi = 2 * PI * ran(0);  /* Scatter randomly in transverse plane */
      cphi = cos(phi);
      sphi = sin(phi);

      rnm = ran(0);	      /* Scatter with H-G dist in prop. dir */

      foo = (1.0 - SQR(gg)) / (1.0 - gg + 2.0 * gg * rnm);
      foo = SQR(foo);
      ctheta = (1.0 + SQR(gg) - foo) / (2.0 * gg);

#if (PARANOID)
      assert((ctheta >= -1) && (ctheta <= 1));
#endif

      theta = acos(ctheta);
      stheta = sin(theta);

      if ((*kz > -1) && (*kz < 1))
	{
	  double bar = sqrt(1 - SQR(kzo));
#if (PARANOID)
	  assert((bar > 0) && (bar <= 1.0));
#endif

	  /* Values of bar on the order 1e-6 causes real problems.
	   * Multiply actual equation by bar and renormalize to 1.0
	   * at the end. */

	  *kx = stheta * (kxo * kzo * cphi - kyo * sphi) + ctheta * kxo * bar;
	  *ky = stheta * (kyo * kzo * cphi + kxo * sphi) + ctheta * kyo * bar;
	  *kz = -stheta * cphi * SQR(bar) + ctheta * kzo * bar;

	  nrm = sqrt(SQR(*kx) + SQR(*ky) + SQR(*kz));

	  *kx /= nrm;
	  *ky /= nrm;
	  *kz /= nrm;
	}
      else
	{
	  /* Vertical direction causes problems, treat as a special case */

	  *kx = stheta * cphi;
	  *ky = stheta * sphi;
	  *kz = ctheta * kzo;
	}
    }

#if (PARANOID)
  /* Sanity-check the output direction cosines */

  assert((*kx >= -1) && (*kx <= 1));
  assert((*ky >= -1) && (*ky <= 1));
  assert((*kz >= -1) && (*kz <= 1));

  assert(SQR(*kx) + SQR(*ky) + SQR(*kz) > 0.9999);
  assert(SQR(*kx) + SQR(*ky) + SQR(*kz) < 1.0001);
#endif

  return;
}

/* ********************************************************************
 *
 * Record a photon in the 2-pt table
 *
 * ********************************************************************/

static void score_photon(const struct Config *syscfg,
			 const TISSUE ***tissueType, const MCCOORD *R,
			 double Ttot, TWOPT_T P2pt, TWOPT_T Q2pt,
			 TWOPT_T *II, TWOPT_T *JJ)
{
  int tindex, VoxN;

  if (syscfg->flags.gen_twopt)
    {
      if (IN_ROI_FULL_C(*R, Ttot, syscfg))
	{
	  tindex = (int)((Ttot - syscfg->Tmin) / syscfg->stepT);

#if (PARANOID)
	  assert(tindex < syscfg->nTstep);
#endif

	  /* Compute index into tables */

	  if (syscfg->flags.matlab_order)
	    {
	      /* Use Matlab ordering so I can dump to disk 
	       * the entire block of memory later 
	       */

	      /* Casting to int is an expensive operation, leave as floats as
	       * long as possible and cast only at the end.  Add a tiny bit
	       * before the cast to avoid the inevitable round-off problems */

	      VoxN = (R->iy - syscfg->Iymin) +
                      syscfg->nIystep * ((R->ix - syscfg->Ixmin) +
                      syscfg->nIxstep * ((R->iz - syscfg->Izmin) +
                      syscfg->nIzstep * (tindex)));
	    }
	  else
	    {
	      /* Use normal C ordering */

	      /* Casting to int is an expensive operation, leave as floats as
	       * long as possible and cast only at the end.  Add a tiny bit 
	       * before the cast to avoid the inevitable round-off problems */

	      VoxN = (R->ix - syscfg->Ixmin) +
                      syscfg->nIxstep * ((R->iy - syscfg->Iymin) +
                      syscfg->nIystep * ((R->iz - syscfg->Izmin) +
                      syscfg->nIzstep * (tindex)));
	    }

#if (PARANOID)
	  assert(VoxN < (syscfg->nIxstep * syscfg->nIystep *
			 syscfg->nIzstep * syscfg->nTstep));
          assert(VoxN >= 0);
#endif

	  if (II != NULL)
	    II[VoxN] += P2pt;

	  if (JJ != NULL)
	    JJ[VoxN] += Q2pt;
	}
    }

  return;
}

/***********************************************************************
 * Advance to the interface.  The interface is on the line segment
 * [x0,y0,z0] -- [x1,y1,z1].  Motion is along the line segment.
 *
 * The algorithm promises to go to an interface.  It does NOT promise
 * to find the _first_ interface.
 **********************************************************************/

static void find_interface(const MCCOORD *_R0, MCCOORD *_R1,
                           const struct Config *system, 
                           const TISSUE ***tissueType)
{
  MCCOORD R, R0, R1;
  double sep;
  TISSUE ti0, ti1, ti;

  /* Make private (modifiable) copies */

  memcpy(&R0, _R0, sizeof(MCCOORD));
  memcpy(&R1, _R1, sizeof(MCCOORD));

  sep = 0;

  /* It doesn't matter what these are, so long as they're different */

  ti0 = TISSUE_TYPE_C(R0, system, tissueType);
  ti1 = TISSUE_TYPE_C(R1, system, tissueType);

  do
    {
      /* Find the midpoint between R0 and R1 */

      assert(ti0 != ti1);

      R.ix = R0.ix;
      R.fx = (R1.ix - R0.ix)/2.0 + (R0.fx + R1.fx) / 2.0;

      R.iy = R0.iy;
      R.fy = (R1.iy - R0.iy)/2.0 + (R0.fy + R1.fy) / 2.0;

      R.iz = R0.iz;
      R.fz = (R1.iz - R0.iz)/2.0 + (R0.fz + R1.fz) / 2.0;

      /* Re-balance fractional portions */

      renorm(&R);

      ti = TISSUE_TYPE_C(R, system, tissueType);

      if (ti == ti0)
        {
          memcpy(&R0, &R, sizeof(MCCOORD));
          ti0 = ti;
        }
      else
        {
          /* Accept, even if ti != ti1 */

          memcpy(&R1, &R, sizeof(MCCOORD));
          ti1 = ti;
        }

      assert(ti0 != ti1);

      sep = DISTANCE(R0,R1);
    }
  while (sep > 1e-6);

  /* Copy back */

  renorm(&R1);
  memcpy(_R1, &R1, sizeof(MCCOORD));

  return;
}

/* Given a fiber NA and a trial photon angle sin(th) (0 is straight
   ahead), accept (returns 1) or reject (returns 0) the angle based on
   a selection from a Gaussian distribution with FWHM of 2*NA. */

static int fiberNA(double NA, double sinth)
{
  double a, b;

#if (HAVE_ASSERT)
  assert((sinth >= -1) && (sinth <= 1));
#endif

  if (NA <= 0.0 || NA > 1)
    return 0;           /* Delta-function distributions, will never match */

  /* Normalization for distribution */

  b = SQR(NA) / LN2;
  a = sqrt(b/PI) * ERF(1/sqrt(b));

  /* Accept if p < f(x) */

  if (ran(0) < exp(-b * SQR(sinth)) / a)
    return 1;
  else
    return 0;
}