commonmc.c

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

/*********************************************************************
 * Photon transport routines, shared between parallel and non-parallel
 *  versions of the code.  Propagate a single photon until it either
 *  exists the system or runs out of time.
 *********************************************************************/

/*
 * This file is part of tMCimg.
 * 
 * tMCimg is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free
 * Software Foundation; either version 2 of the License, or (at your option)
 * any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "config.h"
#include "tMCimg.h"
#include "mccore.h"

#if (HAVE_ASSERT_H)
#include <assert.h>
#endif

#include "commonMC.h"

/*********************************************************************
 Function prototypes
 *********************************************************************/

void run_photon(const struct Config *, int, const TISSUE ***, 
                TWOPT_T *, TWOPT_T *, FILE *, int *);

static void renorm(MCCOORD *);
static void rotate_phase(TWOPT_T *, TWOPT_T *, 
                         const struct Config *, const float *);
static void scatter(double, double *, double *, double *, int *);
static int reflect_photon(const struct Config *, const MCCOORD *, 
                          const MCCOORD *, double *, double *, double *, 
                          int *, const TISSUE ***);
static void find_interface(const MCCOORD *, MCCOORD *,
                           const struct Config *, const TISSUE ***);
static void score_photon(const struct Config *, const TISSUE ***, 
                         const MCCOORD *, double, TWOPT_T, TWOPT_T,
                         TWOPT_T *, TWOPT_T *);

extern void write_photon(const struct Config *, FILE *,
			 double, double, double,
			 double, double, double, double, float *);
static int fiberNA(double, double);

/**************************************************************
  RUN A SINGLE PHOTON THROUGH THE SIMULATION
**************************************************************/

void run_photon(const struct Config *syscfg, int ns,
		const TISSUE ***tissueType, TWOPT_T *II, TWOPT_T *JJ,
		FILE *fp, int *seed)
{
  double  x, y, z;	      /* INITIAL PHOTON POSITION */
  MCCOORD R, Rnew;	      /* CURRENT PHOTON POSITION */
  double  kx, ky, kz;	      /* DIRECTION COSINES */
  double  Lresid, Tresid, Ttot, maxstep, step;
  TWOPT_T P2pt, Q2pt;	      /* PHOTON WEIGHT */
  float   lenTiss[MAXTISSUE]; /* LENGTH SPENT IN EACH TISSUE */
  double  rnm;		      /* RANDOM NUMBER */
  int i,  tissueIndex;
  double  stepT, stepL;

  /* SET THE PHOTON WEIGHT TO 1 AND INITIALIZE PHOTON LENGTH PARAMETERS */

  P2pt = 1.0;
  Q2pt = 0.0;

  Ttot   = 0.0;
  Lresid = 0.0;
  Tresid = 0.0;

  /* Seed RNG, if necessary */

  ran(seed);

  /* Initialize the length in each tissue type */

  for (i = 0; i <= MAXTISSUE; i++)
    lenTiss[i] = 0.0;

  /* Initial source position and direction */

  x = syscfg->srcLoc[ns][0];
  y = syscfg->srcLoc[ns][1];
  z = syscfg->srcLoc[ns][2];

  kx = syscfg->srcDir[ns][0];
  ky = syscfg->srcDir[ns][1];
  kz = syscfg->srcDir[ns][2];

  /* report if direction vector wrong length... */

#if (PARANOID)
  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

  /* Allow for finite source radius */

  if (syscfg->srcRad[ns] > 0)
    {
      double dx, dy, dz, mag;
      double nx, ny, nz, nrm;

    vparallel:
      do
	{
	  /* Generate on a cube then reject based on a sphere */
	  dx = 2 * ran(0) - 1;
	  dy = 2 * ran(0) - 1;
          dz = 2 * ran(0) - 1;
	}
      while (SQR(dx) + SQR(dy) + SQR(dz) > 1.0);

      /* Motion must be in plane perpendicular to initial photon
         direction.  Compute this as the cross-product of my vector
         with the initial k-vector, then rescale back to a unit vector
         to maintain the correct angular distribution */

      nx = ky * dz - kz * dy;
      ny = kz * dx - kx * dz;
      nz = kx * dy - ky * dx;

      if ((nrm = SQR(nx) + SQR(ny) + SQR(nz)) == 0)
        goto vparallel;    /* I know, I know... */
      else
        nrm = sqrt(nrm);   /* Magnitude of the vector */

      /* Turn cross-product into a unit vector to set the angle, and
       * then multiply by a radius selected from a uniform
       * distribution 
       */

      do
	{
	  mag = ran(0);
	}
      while (mag <= 0);

      nx *= mag / nrm;
      ny *= mag / nrm;
      nz *= mag / nrm;

      /* Scale up from a unity-scale vector to the user-specified radius */

      x += nx * (syscfg->srcRad[ns]);
      y += ny * (syscfg->srcRad[ns]);
      z += nz * (syscfg->srcRad[ns]);
    }

  /* Allow for finite source NA */

  if (syscfg->srcNA[ns] > 0)
    {
      double th0, dth;
      double ph0, dph;
      double p1, p2, p3;

      assert(syscfg->srcNA[ns] > 0.0 && syscfg->srcNA[ns] <= 1.0);

      /* The original input direction can be definied as
       * the vector (0,0,1) times a pair of rotation matrices
       * characterized by these two angles */

      th0 = acos(kz);
      ph0 = atan2(ky, kx);

      /* Perturbation d\Omega to the propagation direction */
      /*   NA == n sin(th_0) -> th_0 = asin(NA) (n = 1)    */

      dph = 2 * PI * ran(0);

      do 
        {
          /* Actually sin(theta), convert to sine after we pass
             through fiberNA.  This is not a particularly efficient
             way to select photon angles, but I only have to do it at
             the start of the run. */

          dth = 2*ran(0) - 1;
        }
      while (fiberNA(syscfg->srcNA[ns], dth) != 1);

      dth = asin(dth);

      /* Convert perturbation into rectangular coordinates */

      p1 = sin(dth) * cos(dph);
      p2 = sin(dth) * sin(dph);
      p3 = cos(dth);

      /* Rotate perturbed (0,0,1) vector into its original
       * orientation using the same rotation matrices we would have
       * used for the original (a rotation by th about the y-axis,
       * followed by rotation by ph about the z-axis) 
       */

      kx = cos(th0) * p1 + 0 * p2 + sin(th0) * p3;
      ky = 0 * p1 + 1 * p2 + 0 * p3;
      kz = -sin(th0) * p1 + 0 * p2 + cos(th0) * p3;

      p1 = kx;		  /* Copy back to save a bit of storage space */
      p2 = ky;
      p3 = kz;

      kx = cos(ph0) * p1 + sin(ph0) * p2 + 0 * p3;
      ky = -sin(ph0) * p1 + cos(ph0) * p2 + 0 * p3;
      kz = 0 * p1 + 0 * p2 + 1 * p3;
    }

  /* check rotated length... */

#if (PARANOID)
  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

  /* Go from floating point to "fixed" point system */

  PACKCOORD(R, x, y, z);
  renorm(&R);

  /* Get the initial optical properties */

  tissueIndex = TISSUE_TYPE_C(R, syscfg, tissueType);

#if (PARANOID)
  /* make sure photon doesn't start in air... */
  assert(tissueIndex != 0);
  assert(tissueIndex <= syscfg->Ntissue);
#endif

  /* Loop over scattering events until time exceeds the width of the gate
   * or the photon escapes (leaves whole volume or finds itself in air)
   */

  /* Tresid is the time left until I sample the photon distribution again.
   * The interval between samples is fixed (currently, at compile time).
   * Don't collect before the user-specified starting time. */

  Tresid = MAX(syscfg->Tmin, 0.0) + T_SAMPLE;

  while (Ttot < syscfg->Tmax && IN_TISSUE(tissueIndex) &&
         IN_VOLUME_C(R, syscfg))
    {
      /* Calculate new scattering length ---------------------------- */

      do
	{
	  /* rnm must be strictly positive or the log will blow up and
	   * generate a floating-point error on me 
	   */

	  rnm = ran(0);
	}
      while (rnm <= 0.0);

      /* PROPAGATE THE PHOTON --------------------------------------- */

      /* Lresid is an exponentially distributed random number with a mean
       * of 1.0.  It is also is the dimensionless quantity mu*x 
       *
       * Add to previous value to deal with numerical slop. 
       */

      /* I can probably safely ignore nm-scale differences */

      if (Lresid < L_TOL);
        Lresid = -log(rnm);

      while (Lresid >= L_TOL)
	{
	  assert(IN_VOLUME_C(R, syscfg));

	  /* Compute length of this step ---------------------------- */

	  stepT = Tresid * (syscfg->v0 / syscfg->tn[tissueIndex]);
#if 0
	  stepL = Lresid / syscfg->tmus[tissueIndex];
#else
	  if (syscfg->tmus[tissueIndex] > 0)
	    stepL = Lresid / syscfg->tmus[tissueIndex];
	  else
	    stepL = 1e35;	/* Anything larger than system size works */
#endif

	  assert(stepT > 0);
	  assert(stepL > 0);

	  /* Go with whichever limit is smaller */
	  maxstep = MIN(stepT, stepL);

	  /* If this gets bigger than a voxel, it becomes possible to
	     "tunnel" throuh tissue types, which is bad.  Limiting maxstep
	     reduces these problems. */

	  maxstep = MIN(maxstep, 0.5);	/* Half a voxel steps */
          assert(maxstep > 0);

	  /* Do greedy update of the position until something
	     interesting happens.  Since the calculations involved are
	     minimal, we can use relatively small steps through the voxel */

#if (PARANOID)
	  assert(SQR(kx) + SQR(ky) + SQR(kz) != 0);
#endif

          step = 0;

          do
	    {
              /* L_SAMPLE defined in tMCimg.h, L_SAMPLE=0.5 voxel currently */
              step = MIN(step + L_SAMPLE, maxstep);

              Rnew.ix = R.ix;		/* Copy over interger part */
              Rnew.iy = R.iy;
              Rnew.iz = R.iz;

              Rnew.fx = R.fx + kx * step;	/* update fractional */
              Rnew.fy = R.fy + ky * step;
              Rnew.fz = R.fz + kz * step;

              renorm(&Rnew);		/* Renormalize values */

              if (!IN_VOLUME_C(Rnew, syscfg))
                break;

              if (TISSUE_TYPE_C(Rnew, syscfg, tissueType) != tissueIndex)
                break;
	    }
          while (step < maxstep);

	  /* We already know that this step was too big */
	  maxstep = step;

#if (PARANOID)
          assert((step <= maxstep) && (step > 0));
#endif

	  /* Decide what to do with this step ----------------------- */

	  if (!IN_VOLUME_C(Rnew, syscfg))
	    {
	      /* Stepped out of system, photon is officially gone. */

              memcpy(&R, &Rnew, sizeof(MCCOORD));
	    }
          else 
            {
              /* Getting tissue types is expensive.  Find them once
               * and store the results to be re-used below. */

              int oldIdx = tissueIndex;
              int newIdx = TISSUE_TYPE_C(Rnew, syscfg, tissueType);

#if (PARANOID)
              assert((oldIdx >= 0) && (newIdx >= 0));
#endif

              if (oldIdx == newIdx)
                {
                  /* Trivial case, didn't pass an interface */

                  memcpy(&R, &Rnew, sizeof(MCCOORD));
                }
              else if ((syscfg->tn[oldIdx] == syscfg->tn[newIdx]) &&
                       ((syscfg->flags.mirror == 0) || (newIdx != 0)))
                {
                  /* Mostly trivial case, same index of refraction
                   * before and after stepping, so no possibility for