tmcimg.h

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

/* Generic physical constants - the name tMCimg.h is purely
 *  historic at this point.
 */

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

#ifndef __TMCIMG_H
#define __TMCIMG_H

#include "config.h"

/* Defining PARANOID turns on extra assertions */

#if (PARANOID) && (! HAVE_ASSERT)
# warning "Paranoide checks require assert() - disabling paranoid code"
# undef PARANOID
#endif

/***** Important Constants ********************************************/

#define C_VACUUM  2.99792458e11 	/* Exact, mm/s */

/* This is faster than forever re-computing these values */

#if (HAVE_GSL)
# include "gsl/gsl_math.h"

/* GSL insures that all these definitions exist */
# define SQRT3 M_SQRT3
# define SQRT2 M_SQRT2
# define PI    M_PI
# define LN2   M_LN2

/*
# ifdef SQR
#  undef SQR
# endif

# define SQR(x) gsl_pow_2(x)
*/
#define SQR(x) ((x)*(x))
# define ERF(x) gsl_sf_erf(x)
#else	 				/* !HAVE_GSL */

# ifdef M_PI
#  define PI    M_PI
# else
#  define PI    3.14159265358979323846
# endif

# ifdef M_SQRT2
#  define SQRT2 M_SQRT2
# else
#  define SQRT2 1.41421356237309504880
# endif

# ifdef M_SQRT3
#  define SQRT3 M_SQRT3
# else
#  define SQRT3 1.73205080756887729352
# endif

# ifdef M_LN2
#  define LN2   M_LN2
# else
#  define LN2   0.69314718055994530942
# endif

# ifndef SQR
#  define SQR(x) ((x)*(x))
# endif

# if (HAVE_ERF)
#  define ERF(x) erf(x)
   extern double erf(double);
# else
#  error "erf(x) does not exist"
# endif
#endif

#ifndef ABS
# define ABS(x) (((x) > 0) ? (x) : -(x))
#endif

/* Defaults values */

#define MUS_AIR   -999.0
#define MUA_AIR   -999.0
#define IDX_AIR   -999.0
#define   G_AIR   -999.0

#define FILENMLEN   1024
#define MAXDET       255

/* These #defines only have meaning if we're compiling for parallel
   operation */

/* How often (in seconds) to sample the photon distribution.  Just for
   reference, with 1mm voxels, the voxel crossing time is about 3e-12
   seconds. */

#define T_SAMPLE 10e-12

/* How far (in voxels) to step before re-checking the tissue type.
   Must be less than 1.414; the smaller the better [but, alas, the
   slower as well] */
#define L_SAMPLE 0.5

/* Type of tissueType array.  Change to "short" or "int" if more tissue
 *  types are needed, but that will break compatibility with existing
 *  segmentation files.
 *
 *  NOTE: TISSUE must be a signed type since "-1" is used to indicate
 *  undefined values.
 */

#if (WIDE_TISSUES)
typedef short TISSUE;
# define MAXTISSUE  32767 /* Limited by size of TISSUE */
#else
typedef char TISSUE;
# define MAXTISSUE    127 /* Limited by size of TISSUE */
#endif

/* Using floats instead of doubles cuts the storage and disk space
 * in half, but entails a small sacrifice of accuracy */

#if (SMALL_TWOPT)
typedef float  TWOPT_T;
#else
typedef double TWOPT_T;
#endif

/**** Possibly a faster way to move photons around ********************/

struct _coord {
  int    ix, iy, iz;
  double fx, fy, fz;
};

typedef struct _coord MCCOORD;

/***** Useful macros **************************************************/

/* Given a position, return 1 if the position is inside the defined 
 *  volume, otherwise return 0
 *
 * Casting to integers is an expensive operation on Pentium CPU's,
 *  leave as double as long as I possibly can.
 */

#define IN_VOLUME_C(r, sp) \
    ((r).ix >= 0 && (r).ix < (sp)->nxstep && \
     (r).iy >= 0 && (r).iy < (sp)->nystep && \
     (r).iz >= 0 && (r).iz < (sp)->nzstep)

#define IN_VOLUME(x, y, z, sp) \
    ((x) >= 0 && (x) < (double)(sp)->nxstep && \
     (y) >= 0 && (y) < (double)(sp)->nystep && \
     (z) >= 0 && (z) < (double)(sp)->nzstep)

#define IN_VOLUME_V(r, sp)     IN_VOLUME(r[0],r[1],r[2],sp)

/* Same, but for image volume */

#define IN_ROI_C(R, sp) \
    ((R).ix >= (sp)->Ixmin && (R).ix < (sp)->Ixmax && \
     (R).iy >= (sp)->Iymin && (R).iy < (sp)->Iymax && \
     (R).iz >= (sp)->Izmin && (R).iz < (sp)->Izmax)

#define IN_ROI(x, y, z, sp) \
    ((x) >= (double)(sp)->Ixmin && (x) <= (double)(sp)->Ixmax && \
     (y) >= (double)(sp)->Iymin && (y) <= (double)(sp)->Iymax && \
     (z) >= (double)(sp)->Izmin && (z) <= (double)(sp)->Izmax)

/* Two-point tables care about times as well as position */

#define IN_ROI_FULL_C(R, t, sp) \
    (IN_ROI_C(R,sp) && (t) >= (sp)->Tmin && (t) <= (sp)->Tmax)

#define IN_ROI_FULL(x, y, z, t, sp) \
    (IN_ROI(x,y,z,sp) && (t) >= (sp)->Tmin && (t) <= (sp)->Tmax)

/* Return 1 if the tissue index > 0, which corresponds to propagation
 *  through "tissue" (as opposed to "air")
 */

#define IN_TISSUE(ttidx) ((int)ttidx > 0)

#define INDEX_MATCHED_C(r, sp, tt, ti) \
    ((sp)->tn[(int)TISSUE_TYPE_C(r,sp,tt)] == (sp)->tn[(ti)])

#define INDEX_MATCHED(x, y, z, sp, tt, ti) \
    ((sp)->tn[(int)TISSUE_TYPE(x,y,z,sp,tt)] == (sp)->tn[(ti)])

/* Given a position, return the corresponding tissue type.
 * Note sp is _pointer_ to system struct
 *
 *  Returns -1 for positions outside the defined volume
 */

#define TISSUE_TYPE_C(r, sp, tt) \
   (IN_VOLUME_C(r, sp) ? tt[(r).ix][(r).iy][(r).iz] : (TISSUE)-1)

#define TISSUE_TYPE(x, y, z, sp, tt) \
   (IN_VOLUME(x, y, z, sp) ? tt[(int)(x)][(int)(y)][(int)(z)] : (TISSUE)-1)

/* Short-cut for when caller is POSITIVE that we're in the volume */
#define TISSUE_TYPE_IN_VOL(x, y, z, sp, tt) tt[(int)(x)][(int)(y)][(int)(z)]

/* Backwards compatiblity */
#define GET_TISSUE(x,y,z,sp,tt)  TISSUE_TYPE(x,y,z,sp,tt)
#define GET_TISSUE_V(r,sp,tt)    TISSUE_TYPE(r[0],r[1],r[2],sp,tt)

/***** System Configuration *******************************************/

/* Boolean flags that control program behavior */

struct ProgramFlags
{
  char gen_twopt;		  /* Generate .2pt file? */
  char source_move;		  /* Advance sources to boundary */
  char exact_exit_time;		  /* Store exit time or exit gate */
  char mirror;			  /* Mirrored walls? */
  char matlab_order;		  /* Save in Matlab or C matrix order? */
  char detector_move;		  /* Advance detectors to boundary */
  char exiting_move;		  /* Move exiting photons back to interface */

  char vary_courant;              /* Auto-compute Courant number */
  char fd_method;                 /* Finite difference method to use */
  char bcs;                       /* Finite difference, boundary conditions */
  char box_bcs;                   /* Box BC's */
  char grad_optode;               /* Use gradient optode model */
};

/* Store the entire system configuration.  These values should be
 * treated as read-only once the system configuration routines
 * are finished.
 */

struct Config
{
  int NT;			  /* Total number of photons */
  long seed;			  /* Random number seed */
  double v0;			  /* Speed of light (vacuum) */
  double frequency;		  /* RF Frequency */
  double k0;			  /* 2pi/lambda_rf */
  double dl;			  /* conversion from mm to voxels */

  /* ************************************************************ */

  double optode_depth;            /* FD optode_depth (initial) */
  double dt0;                     /* Default time step (in ps) */
  double courant;                 /* Courant number */

  /* ************************************************************ */

  double xstep, ystep, zstep;	        /* Voxel dimensions */
  int nxstep, Ixmin, Ixmax, nIxstep;	/* Image (ROI) dimensions */
  int nystep, Iymin, Iymax, nIystep;
  int nzstep, Izmin, Izmax, nIzstep;

  double Tmin, Tmax, stepT;	  /* Time stepping for recording */
  int nTstep;

  /* ************************************************************ */

  char segFile[FILENMLEN];	  /* Filename for segmentation file */
  int Ntissue;			  /* Tissue properties */

  double tmus[MAXTISSUE];	  /* Optical properties, by type */
  double tmua[MAXTISSUE];
  double tg[MAXTISSUE];
  double tn[MAXTISSUE];

  /* ************************************************************ */

  int nDets;			  /* # of Detectors */

  double detLoc[MAXDET][3];       /* Position */
  double detDir[MAXDET][3];       /* Orientation (unit vector) */
  double detRad[MAXDET];          /* Radius */
  double  detNA[MAXDET];          /* NA */

  /* ************************************************************ */

  /* tMCimg only supports one source (the first one defined) right
   *   now.  Arrays are defined for future expansion and for code
   *   sharing with other programs.
   */

  int nSrcs;

  double srcLoc[MAXDET][3];
  double srcDir[MAXDET][3];
  double srcRad[MAXDET];
  double  srcNA[MAXDET];

  /* ************************************************************ */

  struct ProgramFlags flags;

#if 0			/* Obsolete flags, don't use */
  double srcRad, srcNA;
  double detRad, detNA;
  double xi, yi, zi;
  double cxi, cyi, czi;

  double Lmin, Lmax, stepL;	  /* min/max photon lenght, computed
				   * from min/max propagation times */
#endif
};

#endif /* __TMCIMG_H */