movesd.c

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

/* ------------------------ MOVE SOURCES/DETECTORS------------------- */

/*
 * 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 "config.h"
#include "tMCimg.h"

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

static void move_srcdir(double *, double *, double *, double, double, double,
                        const struct Config *, const TISSUE ***);
static void find_interface(double *, double *, double *,
                           double *, double *, double *,
                           const struct Config *, const TISSUE ***);

void move_sources(struct Config *system, const TISSUE ***tissueType)
{
  double rx, ry, rz;
  double kx, ky, kz;
  int n;

  for (n = 0; n < system->nSrcs; ++n)
    {
      rx = system->srcLoc[n][0];  /* Starting point */
      ry = system->srcLoc[n][1];
      rz = system->srcLoc[n][2];

      kx = system->srcDir[n][0];
      ky = system->srcDir[n][1];
      kz = system->srcDir[n][2];

      move_srcdir(&rx, &ry, &rz, kx, ky, kz, system, tissueType);
#if (DEBUG)
      printf("Moved source to (%f,%f,%f)\n", rx, ry, rz);
#endif

      system->srcLoc[n][0] = rx;
      system->srcLoc[n][1] = ry;
      system->srcLoc[n][2] = rz;
    }

  return;
}

void move_detectors(struct Config *system, const TISSUE ***tissueType)
{
  double rx, ry, rz;
  double kx, ky, kz;
  int n;

  for (n = 0; n < system->nDets; ++n)
    {
      rx = system->detLoc[n][0];  /* Starting point */
      ry = system->detLoc[n][1];
      rz = system->detLoc[n][2];

      kx = system->detDir[n][0];
      ky = system->detDir[n][1];
      kz = system->detDir[n][2];

      if ((kx == 0) && (ky == 0) && (kz == 0))
        printf("Detector %d not moved, no direction specified\n", n);
      else
        {
          move_srcdir(&rx, &ry, &rz, kx, ky, kz, system, tissueType);
#if (DEBUG)
          printf("Moved detector to (%f,%f,%f)\n", rx, ry, rz);
#endif

          system->detLoc[n][0] = rx;
          system->detLoc[n][1] = ry;
          system->detLoc[n][2] = rz;
        }
    }

  return;
}

/* Common code for sources and detectors, easier and less error-prone
   than maintaing separate functions for each */

static void move_srcdir(double *xold, double *yold, double *zold, 
                        double kx, double ky, double kz,
                        const struct Config *system, 
                        const TISSUE ***tissueType)
{
  TISSUE tissueIndex;
  double x, y, z;
  int loop = 0;

  /* find_interface() finds the interface bounded by the two
     supplied points.  Our starting position is one point [rx,ry,rz], 
     still need to find the second */

  tissueIndex = TISSUE_TYPE(*xold, *yold, *zold, system, tissueType);

  x = *xold;
  y = *yold;
  z = *zold;

  if (IN_TISSUE(tissueIndex))
    {
      /* Too far in, move back */

      while ((loop++ < 1000) && IN_TISSUE(tissueIndex))
        {
          x -= kx;
          y -= ky;
          z -= kz;

          tissueIndex = TISSUE_TYPE(x, y, z, system, tissueType);
        }
    }
  else
    {
      /* Not far enough, move in */

      while ((loop++ < 1000) && !IN_TISSUE(tissueIndex))
        {
          x += kx;
          y += ky;
          z += kz;

          tissueIndex = TISSUE_TYPE(x, y, z, system, tissueType);
        }
    }
      
  if (loop >= 1000)
    {
      fprintf(stderr, "Error finding interface (initial location %f,%f,%f)\n",
              *xold, *yold, *zold);
      exit(1);
    }

  find_interface(xold, yold, zold, &x, &y, &z, system, tissueType);

  /* Sources and detectors should end up INSIDE the tissue */

  if (IN_TISSUE(TISSUE_TYPE(x, y, z, system, tissueType)))
    {
      *xold = x;
      *yold = y;
      *zold = z;
    }

#if (HAVE_ASSERT)
  assert(IN_TISSUE(TISSUE_TYPE(*xold, *yold, *zold, system, tissueType)));
#endif

  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(double *x0, double *y0, double *z0,
                           double *x1, double *y1, double *z1,
                           const struct Config *system, 
                           const TISSUE ***tissueType)
{
  double x, y, z, sep;
  TISSUE ti0, ti1, ti;

  sep = 0;

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

  ti0 = TISSUE_TYPE(*x0, *y0, *z0, system, tissueType);
  ti1 = TISSUE_TYPE(*x1, *y1, *z1, system, tissueType);

#if (HAVE_ASSERT)
  assert(ti0 != ti1);
#endif

  do
    {
      assert(ti0 != ti1);

      x = (*x1 + *x0) / 2;
      y = (*y1 + *y0) / 2;
      z = (*z1 + *z0) / 2;

      ti = TISSUE_TYPE(x, y, z, system, tissueType);

      if (ti == ti0)
        {
          *x0 = x;
          *y0 = y;
          *z0 = z;

          ti0 = TISSUE_TYPE(x, y, z, system, tissueType);
        }
      else
        {
          /* Accept, even if ti != ti1 */

          *x1 = x;
          *y1 = y;
          *z1 = z;

          ti1 = TISSUE_TYPE(x, y, z, system, tissueType);
        }

      assert(ti0 != ti1);

      sep = sqrt(SQR(*x1 - *x0) + SQR(*y1 - *y0) + SQR(*z1 - *z0));
    }
  while (sep > 1e-6);

  return;
}