commonmc.c

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

                   * reflection.  Always accept these steps.
                   *
                   * Checks to avoid special case of mirrored BC's and
                   * the final tissue type of zero--this special case
                   * should fall through into the main reflection code
                   * below.  */

                  memcpy(&R, &Rnew, sizeof(MCCOORD));
                }
              else
                {
                  MCCOORD  Rold;

                  /* Hit a tissue interface, reflect (or not) ----------- */

#if (PARANOID)
                  /* K is supposed to be a unit vector! */
                  assert((SQR(kx)+SQR(ky)+SQR(kz) > 0.9999) &&
                         (SQR(kx)+SQR(ky)+SQR(kz) < 1.00001));

                  assert(oldIdx != newIdx);
                  assert(syscfg->flags.mirror ||
                         (syscfg->tn[oldIdx] != syscfg->tn[newIdx]));
#endif

                  /* Track back to boundary to keep book-keeping straight */

                  memcpy(&Rold, &Rnew, sizeof(MCCOORD));

                  find_interface(&R, &Rnew, syscfg, tissueType);

                  /* Update step to reflect change in location */

                  step -= DISTANCE(Rnew, Rold);

#if (PARANOID)
                  assert(TISSUE_TYPE_C(Rnew, 
                                       syscfg, tissueType) != tissueIndex);

                  assert(IN_VOLUME_C(Rnew, syscfg));
                  /*
                  if (syscfg->flags.mirror)
                    assert(TISSUE_TYPE_C(Rnew, syscfg, tissueType) != 0);
                  */
#endif

                  /* reflect_photon() handles mirrored BC's for me */

                  renorm(&R);
                  renorm(&Rnew);

                  if (reflect_photon(syscfg, &R, &Rnew,
                                     &kx, &ky, &kz, seed, tissueType))
                    {
                      /* Photon reflected, k[xyz] already updated */

#if (PARANOID)
                      assert(TISSUE_TYPE_C(R, syscfg, tissueType)     != 0);
#endif

                      memcpy(&Rold, &Rnew, sizeof(MCCOORD));
                      memcpy(&Rnew, &R,    sizeof(MCCOORD));

                      /* Find the other side of the interface */

                      find_interface(&Rold, &Rnew, syscfg, tissueType);

                      /* Update step to reflect change in location */

                      step -= DISTANCE(Rnew, Rold);

                      /* Now on correct side */
                      memcpy(&R, &Rnew, sizeof(MCCOORD));

#if (PARANOID)
                      assert(TISSUE_TYPE_C(R, syscfg, tissueType) == oldIdx);
#endif
                    }
                  else
                    {
                      /* Photon refracted, k[xyz] already updated ------- */

                      /* Already on correct side */
                      memcpy(&R, &Rnew, sizeof(MCCOORD));
                    }
                }
	    }

#if (PARANOID)
	  assert(IN_TISSUE(tissueIndex));
#endif

	  /* Update based on actual step taken: time ---------------- */

	  Ttot   += step / (syscfg->v0 / syscfg->tn[tissueIndex]);
	  Tresid -= step / (syscfg->v0 / syscfg->tn[tissueIndex]);

	  /* Update based on actual step taken: distance ------------ */

	  Lresid -= step * syscfg->tmus[tissueIndex];

	  lenTiss[tissueIndex] += (float)step;

	  /* If Tresid is zero, score the photon */

	  if (Tresid < T_TOL)
	    {
	      /* record photon in 2pt array, unless out of time */

	      if (Ttot < syscfg->Tmax)
		{
		  /* Compute new amplitude [and phase] from updated lenTiss */
		  rotate_phase(&P2pt, &Q2pt, syscfg, lenTiss);

                  renorm(&R); /* Just in case */

		  score_photon(syscfg, tissueType,
                               &R, Ttot, P2pt, Q2pt, II, JJ);
		}

	      /* Reset Tresid to sample interval, subtract off the
	       * current value just to deal with any minor numerical fuzz.
	       */

	      Tresid = T_SAMPLE - Tresid;
	      assert(Tresid > 0);
	    }

	  /* Update tissue type after making the step */

	  tissueIndex = TISSUE_TYPE_C(R, syscfg, tissueType);

#if (PARANOID)
	  assert(tissueIndex <= syscfg->Ntissue);
#endif

	  /* ----------------------------------------------------------
	   * Break out of the inner loop if we're no longer able to
	   * propagate the photon
	   */

	  if (Ttot < syscfg->Tmax)
            {
#if (DEBUG) && (0)
              printf("Photon out of time\n");
#endif
              break;
            }

	  if (! IN_VOLUME_C(R, syscfg))
            {
#if (DEBUG)
              printf("Photon left imaging volume\n");
#endif
              break;
            }

	  if (! IN_TISSUE(tissueIndex))
            {
#if (DEBUG)
              printf("Photon no longer in tissue\n");
#endif
              break;
            }
	}		      /* Propagate photon */

      /* --------------------------------------------------------------
       * Scatter into new direction using a Heyney-Greenstein
       * distribution if Lresid is zero. 
       */

      if (IN_VOLUME_C(R,syscfg))
        if ((Lresid < L_TOL) && IN_TISSUE(tissueIndex))
          scatter(syscfg->tg[tissueIndex], &kx, &ky, &kz, seed);
    }			      /* Loop until end of single photon */

  /* Score the exiting photon --------------------------------------- */

  if (IN_VOLUME_C(R, syscfg)     && (Ttot < syscfg->Tmax))
    if (!IN_TISSUE(tissueIndex))
      {
	/* Compute new amplitude [and phase] from updated lenTiss */
	rotate_phase(&P2pt, &Q2pt, syscfg, lenTiss);

	/* Save exiting photon to 2pt tables (exit with negative fluence) */

        renorm(&R); /* Just in case */
	score_photon(syscfg, tissueType, &R, Ttot, -P2pt, -Q2pt, II, JJ);

        UNPACKCOORD(R, x, y, z);

        /* Write out photon to .his file.   Parallel and single-threaded
         * versions use different savehis() functions, so call an external
         * wrapper routine here */

        write_photon(syscfg, fp, x, y, z, kx, ky, kz, Ttot, lenTiss);
      }

  return;
}

/* ********************************************************************
 * Tissue changed, have the potential for reflection.
 *
 * Return 1 for reflection, 0 otherwise
 */

static int reflect_photon(const struct Config *syscfg,
                          const MCCOORD *Rold, const MCCOORD *Rnew,
                          double *kx, double *ky, double *kz,
                          int *seed, const TISSUE ***tissueType)
{
  double R_fresnel, R_fresnel_1, R_fresnel_2;
  double ctheta, ctheta_out, ctheta_c2;
  double sn1, sn2, sn3;
  double nold, nnew, rnm;
  int    norm, tissueIndex_old, tissueIndex_new;

#if (PARANOID)
  assert(IN_VOLUME_C(*Rnew, syscfg));
#endif

  tissueIndex_old = TISSUE_TYPE_C(*Rold, syscfg, tissueType);
  tissueIndex_new = TISSUE_TYPE_C(*Rnew, syscfg, tissueType);

  assert(tissueIndex_old != tissueIndex_new);

#if (PARANOID)
  assert(tissueIndex_old >= 0 && tissueIndex_old <= syscfg->Ntissue);
  assert(tissueIndex_new >= 0 && tissueIndex_new <= syscfg->Ntissue);
#endif

  nold = syscfg->tn[tissueIndex_old];
  nnew = syscfg->tn[tissueIndex_new];

  /* Change place-holder values into 1.0 */

  if (nold == IDX_AIR)
    nold = 1.0;
  if (nnew == IDX_AIR)
    nnew = 1.0;

  /* 2.75 is physically reasonable, but arbitrary */
#if (PARANOID)
  assert((nold >= 1.0) && (nold < 2.75));
  assert((nnew >= 1.0) && (nnew < 2.75));
#endif

  /* CALCULATE SURFACE NORMAL VECTOR OUT OF TISSUE */

  /* There are really only 3 possible cases using
   * our linearized approach
   *
   * 1) Along a single coordinate axis (cube face)
   * 2) Pair of coordinate axis (cube edge)
   * 3) Along 3 coordinate axes (cube corner)
   * 
   * Take the interface to be orthogonal to one
   * of these three instances 
   */

  /* THIS IS SERIOUSLY BROKEN! REPORTS CORNERS IF WE CROSS DIAGONALLY
     INTO AN ADJACENT VOXEL, EVEN THOUGH THE INTERFACE IS PLANAR! */

  sn1 = Rnew->ix - Rold->ix;
  sn2 = Rnew->iy - Rold->iy;
  sn3 = Rnew->iz - Rold->iz;

  /* inew-i = 0 or +/- 1, so the square is always 0 or 1 */
  norm = SQR(Rnew->ix - Rold->ix) +  
    SQR(Rnew->iy - Rold->iy) + SQR(Rnew->iz - Rold->iz);

  /* Because its computed from integers, the normalization
   * factors can be specified excatly in advance - saves time */

  switch (norm)
    {
      case 1:
	break;			  /* Already normalized */
      case 2:
	sn1 /= SQRT2;
	sn2 /= SQRT2;
	sn3 /= SQRT2;
	break;
      case 3:
	sn1 /= SQRT3;
	sn2 /= SQRT3;
	sn3 /= SQRT3;
	break;
      default:
	fprintf(stderr, "Program bug, norm = %d != 1,2, or 3\n", norm);
#if (PARANOID)
	assert((norm == 1) || (norm == 2) || (norm == 3));
#endif
	exit(1);
    }

  /* cosine of angle of incidence */
  ctheta = (*kx) * sn1 + (*ky) * sn2 + (*kz) * sn3;

  /* squared cosine of critical angle (incident) */
  ctheta_c2 = 1 - SQR(nnew) / SQR(nold);

  if (SQR(ctheta) < ctheta_c2)
    {
      /* Total internal reflection, always reflected */

      R_fresnel = 1;
      ctheta_out = -999;	/* dummy value */
    }
  else if (syscfg->flags.mirror &&
	   (TISSUE_TYPE_C(*Rnew, syscfg, tissueType) == 0))
    {
      /* Perfectly reflecting air-tissue bound's, always reflected */

      R_fresnel = 1;
      ctheta_out = -999;	/* dummy value */
    }
  else
    {
      /* cosine of transmitted angle */

      ctheta_out = sqrt(1 - (SQR(nold) / SQR(nnew)) * (1 - SQR(ctheta)));

      /* Calculate Fresnel reflection coefficient */

      R_fresnel_1 = ((nold * ctheta_out - nnew * ctheta) /
		     (nold * ctheta_out + nnew * ctheta));
      R_fresnel_2 = ((nold * ctheta - nnew * ctheta_out) /
		     (nold * ctheta + nnew * ctheta_out));

      R_fresnel = (SQR(R_fresnel_1) + SQR(R_fresnel_2)) / 2;
    }

  /* Reflect if ran() < R_fresnel */

  if ((rnm = ran(0)) <= R_fresnel)
    {
      /* Reflection occurred, adjust wave-vector, k.n = -k.n */

#if (PARANOID)
      assert(TISSUE_TYPE_C(*Rold, syscfg, tissueType) != 0);
#endif

      *kx -= 2 * ctheta * sn1;
      *ky -= 2 * ctheta * sn2;
      *kz -= 2 * ctheta * sn3;

#if (PARANOID)
      assert((*kx >= -1) && (*kx <= 1));
      assert((*ky >= -1) && (*ky <= 1));
      assert((*kz >= -1) && (*kz <= 1));
      assert(SQR(*kx) + SQR(*ky) + SQR(*kz) > 0.9999);
      assert(SQR(*kx) + SQR(*ky) + SQR(*kz) < 1.0001);
#endif

      return 1;
    }
  else
    {
      double knorm;

      /* Transmission occured, refract according to Snell's law */

#if (PARANOID)
      assert((ctheta_out >= -1) && (ctheta_out <= 1));
#endif

      /* Subtract cos(theta), leaving sin(theta) terms */
      *kx -= ctheta * sn1;
      *ky -= ctheta * sn2;
      *kz -= ctheta * sn3;

      /* Rescale sin(theta), preserving sin(phi),cos(phi)
       * to get sin(theta') [E_{||} terms] */

      *kx *= nold / nnew;
      *ky *= nold / nnew;
      *kz *= nold / nnew;

      /* Add back in cos(theta) [E_\perp terms] */

      *kx += ctheta_out * sn1;
      *ky += ctheta_out * sn2;
      *kz += ctheta_out * sn3;

#if (PARANOID)
      assert((*kx >= -1) && (*kx <= 1));
      assert((*ky >= -1) && (*ky <= 1));
      assert((*kz >= -1) && (*kz <= 1));
#endif

      /* Re-normalize to keep K a unit vector */

      knorm = sqrt(SQR(*kx) + SQR(*ky) + SQR(*kz));
      *kx /= knorm;