pmcimg.c

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

/*********************************************************************
 * pMCimg
 *
 * Parallelized version of tMCimg.
 *
 *********************************************************************/

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

#if (HAVE_UNISTD_H)
#include <unistd.h>
#endif

#if (! HAVE_PERROR)
/* Easier than endless #ifdef's */
# define perror(x) fprintf(stderr, "%s\n", x)
#endif

#ifndef USE_MPI
# error "MPI not requested, but trying to build parallel code"
#else
/* Get defines specific to the parallel code */
# include "parallel.h"
#endif

#define USAGE printf("optionss: 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 source 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 *, 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 *);

extern char *mpibuff;

/*********************************************************************
 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 */
  TWOPT_T *II0,*JJ0;	      /* Extra copy for master program */
  TISSUE ***tissueType;	      /* STORE THE IMAGE FILE */
  int i, j, n, nvox;

  /* INITIALIZE THE CONFIGURATION STRUCTURE'S MEMORY */

  initialize_config(&syscfg);

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

  fpusetup();

  /* SET UP THE MPI LIBRARY */

  initParallel(&argc, &argv);	/* Call before parsing arguments! */

#if 0
  if (getNodeCount() <= 2)
    {
      fprintf(stderr, 
              "Parallelized code does not support nodecount %d <= 2\n",
	      getNodeCount());
      fprintf(stderr, "Call tMCimg for single-threaded simulations instead\n");
      MPI_Abort(MPI_COMM_WORLD, 1);
    }
#endif

  /* 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)
    MPI_Abort(MPI_COMM_WORLD, n);

  if (syscfg.NT < 1)
    {
      fprintf(stderr, "Error: total number of photons must be >1 [is %d]\n",
	syscfg.NT);
      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)
    {
      MPI_Abort(MPI_COMM_WORLD, n);
      exit(n);			  /* probably not reached */
    }

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

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

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

      if ((II0 = alloc_1d_array(nvox,sizeof(TWOPT_T))) == NULL)
	MPI_Abort(MPI_COMM_WORLD, 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)
	    MPI_Abort(MPI_COMM_WORLD, 1);
	  if ((JJ0 = alloc_1d_array(nvox,sizeof(TWOPT_T))) == NULL)
	    MPI_Abort(MPI_COMM_WORLD, 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;

      if (II != NULL && II0 != NULL)
	memcpy(II0, II, nvox*sizeof(TWOPT_T));

      if (JJ != NULL && JJ0 != NULL)
	memcpy(JJ0, JJ, nvox*sizeof(TWOPT_T));
    }