hisfile.c

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

/********************************************************************
 SAVE PHOTON TO HISTORY FILE : this converts length back to mm units.
*********************************************************************/

/*
 * 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_FCNTL_H)
#include <fcntl.h>
#endif

#if (! HAVE_PERROR)
/* Use fprintf instead of perror.  Not as nice, but it will suffice. */
#define perror(s)  fprintf(stderr,s)
#endif

/*
extern float unscale_length(const float, const struct Config *);
*/

/* OPEN A FILE POINTER TO SAVE THE HISTORY INFORMATION */

FILE *openhis(int ns, const char *basefile, const struct Config *system)
{ 
  char filenm[FILENMLEN];
  FILE *fp;

#if (HAVE_SNPRINTF)
  if (system->nDets == 0)
    {
      if (system->nSrcs > 1)
        snprintf(filenm, FILENMLEN, "%s.his.a%d", basefile, ns);
      else
        snprintf(filenm, FILENMLEN, "%s.his.a", basefile);
    } 
  else if (system->flags.exact_exit_time != 0) 
    {
      if (system->nSrcs > 1)
        snprintf(filenm, FILENMLEN, "%s.his.t%d", basefile, ns);
      else
        snprintf(filenm, FILENMLEN, "%s.his.t", basefile);
    } 
  else 
    {
      if (system->nSrcs > 1)
        sprintf(filenm, "%s.his.%d", basefile, ns);
      else
        sprintf(filenm, "%s.his", basefile);
    }
#else  /* HAVE_SNPRINTF */  
  if (system->nDets == 0)
    {
      if (system->nSrcs > 1)
        sprintf(filenm, "%s.his.a%d", basefile, ns);
      else
        sprintf(filenm, "%s.his.a", basefile);
    } 
  else if (system->flags.exact_exit_time != 0) 
    {
      if (system->nSrcs > 1)
        sprintf(filenm, "%s.his.t%d", basefile, ns);
      else
        sprintf(filenm, "%s.his.t", basefile);
    } 
  else 
    {
      if (system->nSrcs > 1)
        sprintf(filenm, "%s.his.%d", basefile, ns);
      else
        sprintf(filenm, "%s.his", basefile);
    }
#endif

  /* Write or append?  Write is correct for single-CPU implementations,
   *  parallel implementations require append.
   */

  if ((fp = fopen(filenm, "wb")) == NULL) 
    {
      fprintf(stderr, "ERROR: Can't open history file for writing\n");
      exit(1);
    }

  return fp;
}

void savehis(const struct Config *system, FILE *fp,
	     double x, double y, double z,
	     double kx, double ky, double kz, double T, float *lenTiss)
{
  float buff[4];		  /* temporary storage */
  int i, tindex;
  size_t nf;

  /* Current time gate */

  tindex = (int)((T - system->Tmin) / system->stepT);

  /* Rescale back to physical units */

  for (i = 1; i <= system->Ntissue; i++)
    lenTiss[i] = unscale_length(lenTiss[i], system);

  /* IF NO DETECTORS DEFINED THEN SAVE EXIT POSITION 
   * ELSE LOOP THROUGH NUMBER OF DETECTORS */

  if (system->nDets == 0)
    {
      /* Convert back to normal dimensions [mm] before saving */

      buff[0] = unscale_length((float)x, system);
      buff[1] = unscale_length((float)y, system);
      buff[2] = unscale_length((float)z, system);

      if (system->flags.exact_exit_time != 0)
	buff[3] = (float)T;
      else
	buff[3] = (float)tindex;

      if ((nf = fwrite(&buff, sizeof(float), 4, fp)) != 4)
	{
	  perror("Error writing to history file");
	  exit(1);
	}

      if ((nf = fwrite(&lenTiss[1], 
                       sizeof(float), (size_t)system->Ntissue, fp)) 
          != (size_t)system->Ntissue)
	{
	  perror("Error writing to history file");
	  exit(1);
	}
    }
  else				  /* system->nDets > 0 */
    {
      for (i = 0; i < system->nDets; i++)
	{
	  /* dR and system->detRad[] are both in voxel units */

	  double dR = SQR(x - system->detLoc[i][0]) +
	              SQR(y - system->detLoc[i][1]) + 
                      SQR(z - system->detLoc[i][2]);

	  if (dR <= SQR(system->detRad[i]))
	    {
	      buff[0] = (float) i;	/* Detector #  */

	      /* (kx,ky,kz) is implicitly a unit vector */

#if (HAVE_ASSERT)
	      assert((SQR(kx) + SQR(ky) + SQR(kz) > 0.9999) &&
		     (SQR(kx) + SQR(ky) + SQR(kz) < 1.0001));
#endif

	      /* Accept/reject based on detector NA */
	      /* Accept both NA=0 and NA=1 to specify all photons
	       *  should be accepted (i.e. effectively disable the check) */

	      if ((system->detNA[i] > 0) && (system->detNA[i] < 1))
		{
		  double p1, p2, p3, proj;

		  /* normals are implicitly unit vectors */

		  printf("system->detNA[%d]=%f\n", i, system->detNA[i]);
#if (HAVE_ASSERT)
		  assert(SQR(system->detDir[i][0]) +
			 SQR(system->detDir[i][1]) +
			 SQR(system->detDir[i][2]) > 0.9999);
		  assert(SQR(system->detDir[i][0]) +
			 SQR(system->detDir[i][1]) +
			 SQR(system->detDir[i][2]) < 1.0001);
#endif

		  /* Resolve the sign ambiguity in the normals (due
		   * to the different possible sign conventions) by
		   * only working with the magnitude of the k.n.
		   *
		   * Photons moving back into tissue from air aren't
		   *  allowed, so this should be an acceptable solution.
		   */

		  p1 = kx * system->detDir[i][0];
		  p2 = ky * system->detDir[i][1];
		  p3 = kz * system->detDir[i][2];

		  /* Find sin() of angle between detDir and K, computed
		   *   as sqrt(1 - cos()^2) where the cosine is also the
		   *   dot product */
		  proj = sqrt(1 - (SQR(p1) + SQR(p2) + SQR(p3)));

		  /* Compare projected angle to NA of detector */
		  if (proj > system->detNA[i])
		    continue;	  /* Out of acceptance angle */
		}

	      /* use either tgate # or exact photon exit time... */

	      if (system->flags.exact_exit_time != 0)
		buff[1] = (float) T;
	      else
		buff[1] = (float) tindex;

	      /* WRITE TO THE HISTORY FILE */

	      if ((nf = fwrite(&buff, sizeof(float), 2, fp)) != 2)
		{
		  perror("Error writing to history file");
		  exit(1);
		}

	      if ((nf = fwrite(&lenTiss[1], sizeof(float),
			       (size_t) system->Ntissue,
			       fp)) != (size_t) system->Ntissue)
		{
		  perror("Error writing to history file");
		  exit(1);
		}
	    }
	}
    }

  return;
}