tmcimg.c

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

/*********************************************************************
 * tMCimg
 *
 *  mods from tMCimg6: barnett 3/1/02 added cmd line options -nst
 *                     fixed voxel<->mm dimensional problems.
 *      HIS formats    .his variants are .his.t (det #, exact times) and
 *                     .his.a (all photon locations, exact times) (nDet=0)
 *  4/4/02 Jon's mirrored-air-bdry mod cmd line, tissuetype update bugfix.
 *  
 *  11Apr02 Merged JS and AB versions, changed filename from
 *            tMCimg6.c to tMCimg6b.c.  Moved flags into the config
 *            structure.  Changed parser to allow different radii, NA,
 *            direction, and position for every source and detector.
 *            tMCimg still only uses the first source and does not yet
 *            support detector NA though.
 *
 *  March 2003 Created tMCimg8.c, which solves the previously
 *            identified problem with sampling unevenly near 
 *            interfaces.  Lots of other small fixes.
 *
 *********************************************************************/

/*
 * 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"

#define USAGE printf("options: if string contains...\n"\
          "\t? - display the option help (this)\n"\
          "\tn - inhibit .2pt file storage and dump\n"\
          "\tR - Use normal C ordering for segmentation/.2pt files\n"\
          "\ts - don't move sources to the boundary\n"\
          "\td - don't move detectors to the boundary\n"\
          "\tt - write exact exit time to .his.t instead of tgate # to .his\n"\
          "\tm - make tissue-air boundaries perfectly mirrored\n\n")

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

static void run_simulation(const struct Config *, const TISSUE ***,
                           TWOPT_T *, TWOPT_T *, const char *);

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

void write_photon(const struct Config *, FILE *, double, double, double,
		  double, double, double, double, float *);

/*********************************************************************
 BEGIN PROGRAM CODE
 *********************************************************************/

int main(int argc, char *argv[])
{
  struct Config syscfg;	      /* System configuration, constructed
			       *   from the .cfg file */
  TWOPT_T *II, *JJ;	      /* FOR STORING THE 2-PT FLUENCE */
  TISSUE ***tissueType;	      /* STORE THE IMAGE FILE */
  int i, j, n, nvox;

  /* Initialize the configuration structure's memory */

  initialize_config(&syscfg);

  /* GET THE COMMAND LINE ARGUMENT(S) */

  if (argc == 1)
    {
      printf("tMCimg version 8.0%c\n", 'A' + 0);
      printf("USAGE:: tMCimg [-opts] input_file (.cfg assumed)\n");

      USAGE;

      printf("Default options:\n");
      printf("\t2pt files %s be generated\n",
	     (syscfg.flags.gen_twopt != 0) ? "will" : "will not");
      printf("\t%s-ordered data files\n",
	     (syscfg.flags.matlab_order != 0) ? "Matlab" : "C");
      printf("\tSources %s to interface\n",
	     (syscfg.flags.source_move != 0) ? "moved" : "not moved");
      printf("\tDetectors %s to interface\n",
	     (syscfg.flags.detector_move != 0) ? "moved" : "not moved");
      printf("\tExiting photons %s to interface\n",
	     (syscfg.flags.exiting_move != 0) ? "moved" : "not moved");
      printf("\t%s recorded in history file\n",
	     (syscfg.flags.exact_exit_time != 0) ?
	     "Exact exit times" : "Time indices");
      if (syscfg.flags.mirror != 0)
	printf("\tMirrored air boundaries\n");
      printf("\n");

      return 1;
    }
  else
    for (i = 1; i < argc; i++)
      {
	if (argv[i][0] == '-')	/* introduces option flags */
	  {
	    /* search through command line option words to find
	     * meaningful chars... */

	    for (j = 1; argv[i][j] != '\0'; j++)
	      switch (argv[i][j])
		{
                case 'm':
                  syscfg.flags.mirror = !syscfg.flags.mirror;
                  printf("Tissue-air boundaries %s be mirrored.\n",
                         (syscfg.flags.mirror ? "will" : "will not"));
                  break;
                case 'n':
                  syscfg.flags.gen_twopt = !syscfg.flags.gen_twopt;
                  printf("The .2pt file %s be generated.\n",
			 (syscfg.flags.gen_twopt ? "will" : "will not"));
                  break;
                case 'R':
                  syscfg.flags.matlab_order = !syscfg.flags.matlab_order;
                  printf("Data %s be saved in matlab order.\n",
                         (syscfg.flags.matlab_order ? "will" : "will not"));
                  break;
                case 'd':
                  syscfg.flags.detector_move = !syscfg.flags.detector_move;
                  printf("Detectors %s be moved to interface.\n",
                         (syscfg.flags.detector_move ? "will" : "will not"));
                  break;
                case 's':
                  syscfg.flags.source_move = !syscfg.flags.source_move;
                  printf("Source %s be moved to interface.\n",
                         (syscfg.flags.source_move ? "will" : "will not"));
                  break;
                case 'e':
                  syscfg.flags.exiting_move = !syscfg.flags.exiting_move;
                  printf("Exiting photons %s be moved to boundary.\n",
                         (syscfg.flags.exiting_move ? "will" : "will not"));
                  break;
                case 't':
                  syscfg.flags.exact_exit_time = !syscfg.flags.exact_exit_time;
                  printf("History file %s record exact exit times.\n",
                         (syscfg.flags.exact_exit_time ? "will" : "will not"));
                  break;
                case 'h':
                case '?':
                  USAGE;
                  return 0;
                default:
                  printf("Unknown option -%c ignored\n", argv[i][j]);
		}
	  }
	else
	  {
	    /* Assume that since it isn't a command-line flag that it's
	     *  the name of the configuration file.  Diddle the argument
	     *  vector to fool the parser into doing the right thing.
	     */

	    argv[1] = argv[i];
	    argv[2] = NULL;   /* may or may not be necessary */

	    argc = 2;

	    break;
	  }
      }

  /* PARSE THE INPUT FILE */

  if ((n = loadconfig(&syscfg, argv[1])) != 0)
    exit(n);

  if (syscfg.NT < 1)
    {
      fprintf(stderr, "Error: total number of photons must be >1\n");
      exit(1);
    }

  rescale_system(&syscfg);

  /* SUMMARY OF EXTANT OPTIONS */

  if (syscfg.flags.gen_twopt)
    printf("...2pt files will be generated\n");
  if (syscfg.flags.matlab_order)
    printf("...Matlab-ordered data files\n");
  if (syscfg.flags.source_move)
    printf("...Sources moved to interface\n");
  if (syscfg.flags.detector_move)
    printf("...Detectors moved to interface\n");
  if (syscfg.flags.exiting_move)
    printf("...Exiting photons moved to interface\n");
  if (syscfg.flags.exact_exit_time)
    printf("...Exact exit times recorded in history file\n");
  if (syscfg.flags.mirror)
    printf("...Mirrored air boundaries\n");

  /* READ IN THE SEGMENTED DATA FILE */

  if ((n = loadsegments(&syscfg, &tissueType)) != 0)
    exit(n);

  /* ALLOCATE SPACE FOR AND INITIALIZE THE PHOTON FLUENCE TO 0 */

  nvox = syscfg.nIxstep * syscfg.nIystep * syscfg.nIzstep * syscfg.nTstep;
  II = NULL;
  JJ = NULL;

  if (syscfg.flags.gen_twopt)
    {
      if ((II = alloc_1d_array(nvox,sizeof(TWOPT_T))) == NULL)
	exit(1);

      /* CW/TD don't need this.  Hopefully, the viewer etc. will be
       *  smart enought to know the difference */

      if (syscfg.frequency != 0.0)
	if ((JJ = alloc_1d_array(nvox,sizeof(TWOPT_T))) == NULL)
	  exit(1);

      /* INITIALIZE THE FLUENCE DATA */

      if (II != NULL)
	for (i = 0; i < nvox; i++)
	  II[i] = 0.0;

      if (JJ != NULL)
	for (i = 0; i < nvox; i++)
	  JJ[i] = 0.0;
    }

  /* SET UP THE FLOATING POINT CO-PROCESSOR */

  fpusetup();

  /* MAKE SURE THE SOURCES ARE AT INTERFACES */

  if (syscfg.flags.source_move)
    move_sources(&syscfg, (const TISSUE ***)tissueType);

  /* MAKE SURE THE DETECTORS ARE AT INTERFACES */

  if (syscfg.flags.detector_move)
    move_detectors(&syscfg, (const TISSUE ***)tissueType);

  /* RUN THE MONTE CARLO SIMULATION */

  run_simulation(&syscfg, (const TISSUE ***)tissueType, II, JJ, argv[1]);

  /* DONE */

  nvox = syscfg.nIxstep * syscfg.nIystep * syscfg.nIzstep * syscfg.nTstep;

  free_1d_array(II, nvox);
  free_1d_array(JJ, nvox);

  return 0;
}

/* ------------------------ RUN THE SIMULATION ---------------------- */

static void run_simulation(const struct Config *syscfg, 
                           const TISSUE ***tissueType,
                           TWOPT_T *II, TWOPT_T *JJ, const char *basefile)
{
  FILE *fp = NULL;		  /* FILE POINTERS FOR SAVING THE DATA */
  int i, ns, nvox, seed;
  unsigned long N;
#if (REALLY_RANDOM)
  FILE *rp;
#endif

#if (REALLY_RANDOM)
  /* /dev/random is a much better seed for parallel code than seed+n */

  if ((rp = fopen("/dev/random", "rb")) == NULL)
    {
      if ((rp = fopen("/dev/urandom", "rb")) == NULL)
	{
	  fread(&seed, sizeof(seed), 1, rp);
	  fclose(rp);
	}
      else
	{
	  printf("Warning: can't open /dev/random, using less random seed\n");
	  seed += myNode();
	}
    }
  else
    {
      fread(&seed, sizeof(seed), 1, rp);
      fclose(rp);
    }
#else
  seed = syscfg->seed;
#endif

  nvox = syscfg->nIystep * syscfg->nIxstep * syscfg->nIzstep * syscfg->nTstep;

  if (syscfg->nSrcs > 1)
    {
      if (syscfg->nDets == 0)
	printf("No detectors: writing all exiting photons to %s.his.aN\n",
	       basefile);
      else if (syscfg->flags.exact_exit_time)
	printf("Writing exact exit times to %s.his.tN\n", basefile);
      else
	printf("Writing exiting tgate number to %s.his.N\n", basefile);
    }
  else
    {
      if (syscfg->nDets == 0)
	printf("No detectors: writing all exiting photons to %s.his.a\n",
	       basefile);
      else if (syscfg->flags.exact_exit_time)
	printf("Writing exact exit times to %s.his.t\n", basefile);
      else
	printf("Writing exiting tgate number to %s.his\n", basefile);
    }


  /* FOR EVERY SOURCE, RUN THE SIMULATION */

  for (ns = 0; ns < syscfg->nSrcs; ns++)
    {
      fp = openhis(ns, basefile, syscfg);

      if (fp == NULL)
        {
          fprintf(stderr, "Can't open history file %s\n", basefile);
          exit(1);
        }

      /*******************************************************************
       * START MIGRATING THE PHOTON GENERATING PHOTONS UNTIL NUMBER OF
       * PHOTONS EXECUTED (N) IS EQUAL TO THE NUMBER TO BE GENERATED (NT)
       *******************************************************************/

      for (N = 0; N < (unsigned long)syscfg->NT; N++)
        {
	  if (syscfg->NT >= 10 && ((N+1) % (syscfg->NT / 10)) == 0)
	    {
	      printf("Completed %ld of %d photons\n", N+1, syscfg->NT);
	      fflush(stdout);
	    }

	  run_photon(syscfg, ns, tissueType, II, JJ, fp, &seed);
	}

      /* CLOSE HISTORY FILE */

      fclose(fp);

      /* SAVE THE PHOTON DENSITY (2-POINT GREENS FUNCTION) */

      if (syscfg->flags.gen_twopt)
	{
	  /* Print some useful information if you're manipulating tables
	   *  by hand from inside matlab */

	  printf("2pt file = %d x %d x %d x %d\n",
		 syscfg->nIystep, syscfg->nIxstep,
		 syscfg->nIzstep, syscfg->nTstep);

	  /* Save fluences to disk */

	  if (syscfg->nSrcs == 1)
	    save2pt(basefile, -1, II, JJ, nvox, NULL);
	  else
	    {
	      save2pt(basefile, ns, II, JJ, nvox, NULL);

	      /* Reset memory for next pass through */

	      if (II != NULL)
		for (i = 0; i < nvox; i++)
		  II[i] = 0.0;

	      if (JJ != NULL)
		for (i = 0; i < nvox; i++)
		  JJ[i] = 0.0;
	    }
	}
    }

  return;
}

/* Call the single-threaded savehis() routine */

void write_photon(const struct Config *syscfg, FILE *fp,
		  double x, double y, double z,
		  double kx, double ky, double kz, 
		  double Ttot, float *lenTiss)
{
  savehis(syscfg, fp, x, y, z, kx, ky, kz, Ttot, lenTiss);

  return;
}