Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
MitkMCxyz.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 // Please retain the following copyright notice
13 /******************************************************************
14  * based on mcxyz.c Oct2014
15  *
16  * mcxyz.c, in ANSI Standard C programing language
17  *
18  * created 2010, 2012 by
19  * Steven L. JACQUES
20  * Ting LI
21  * Oregon Health & Science University
22  *
23  *******************************************************************/
24 
25 #include <math.h>
26 #include <stdio.h>
27 #include <string.h>
28 #include <sstream>
29 #include <stdlib.h>
30 #include <time.h>
31 #include <thread>
32 #include <chrono>
33 
34 #include <vector>
35 #include <iostream>
36 
37 #include "mitkCommandLineParser.h"
38 #include "mitkIOUtil.h"
39 #include "mitkImageCast.h"
40 #include <mitkImage.h>
41 #include <mitkImageReadAccessor.h>
42 #include <mitkCoreServices.h>
45 
46 #include <mitkPAProbe.h>
47 #include <mitkPALightSource.h>
49 
50 #ifdef _WIN32
51 #include <direct.h>
52 #else
53 #include <sys/types.h>
54 #include <sys/stat.h>
55 #endif
56 
57 #define ls 1.0E-7 /* Moving photon a little bit off the voxel face */
58 #define PI 3.1415926
59 #define ALIVE 1 /* if photon not yet terminated */
60 #define DEAD 0 /* if photon is to be terminated */
61 #define THRESHOLD 0.01 /* used in roulette */
62 #define CHANCE 0.1 /* used in roulette */
63 #define SQR(x) (x*x)
64 #define SIGN(x) ((x)>=0 ? 1:-1)
65 #define ONE_MINUS_COSZERO 1.0E-12 /* If 1-cos(theta) <= ONE_MINUS_COSZERO, fabs(theta) <= 1e-6 rad. */
66  /* If 1+cos(theta) <= ONE_MINUS_COSZERO, fabs(PI-theta) <= 1e-6 rad. */
67 
68  /* Struct for storing x,y and z coordinates */
69 struct Location
70 {
71  int x, y, z;
72  double absorb;
73 };
74 
75 struct Location initLocation(int x, int y, int z, double absorb)
76 {
77  struct Location loc;
78  loc.x = x;
79  loc.y = y;
80  loc.z = z;
81  loc.absorb = absorb;
82  return loc;
83 }
84 
85 class DetectorVoxel
86 {
87 public:
88  Location location;
89  std::vector<Location>* recordedPhotonRoute = new std::vector<Location>();
90  double* fluenceContribution;
91  double m_PhotonNormalizationValue;
92  long m_NumberPhotonsCurrent;
93 
94  DetectorVoxel(Location location, long totalNumberOfVoxels, double photonNormalizationValue)
95  {
96  this->location = location;
97  this->fluenceContribution = (double *)malloc(totalNumberOfVoxels * sizeof(double));
98  for (int j = 0; j < totalNumberOfVoxels; j++) fluenceContribution[j] = 0; // ensure fluenceContribution[] starts empty.
99  m_NumberPhotonsCurrent = 0;
100  m_PhotonNormalizationValue = photonNormalizationValue;
101  }
102 };
103 
104 bool verbose(false);
105 
106 class InputValues
107 {
108 private:
109  std::string inputFilename;
110  int tissueIterator;
111  long long ix, iy, iz;
112 public:
113  int mcflag, launchflag, boundaryflag;
114  double xfocus, yfocus, zfocus;
115  double ux0, uy0, uz0;
116  double radius;
117  double waist;
118  double xs, ys, zs; /* launch position */
119  int Nx, Ny, Nz, numberOfTissueTypes; /* # of bins */
120 
121  char* tissueType;
122  double* muaVector;
123  double* musVector;
124  double* gVector;
125  double* normalizationVector;
126 
127  double xSpacing, ySpacing, zSpacing;
128  double simulationTimeFromFile;
129  long long Nphotons;
130  long totalNumberOfVoxels;
131  double* totalFluence;
132  std::string myname;
133  DetectorVoxel* detectorVoxel;
134  mitk::Image::Pointer m_inputImage;
135  mitk::Image::Pointer m_normalizationImage;
136 
137  InputValues()
138  {
139  tissueType = nullptr;
140  muaVector = nullptr;
141  musVector = nullptr;
142  gVector = nullptr;
143  detectorVoxel = nullptr;
144  normalizationVector = nullptr;
145 
146  mcflag = 0;
147  launchflag = 0;
148  boundaryflag = 0;
149  }
150 
151  double GetNormalizationValue(int x, int y, int z)
152  {
153  if (normalizationVector)
154  return normalizationVector[z*Ny*Nx + x*Ny + y];
155  else
156  return 1;
157  }
158 
159  void LoadValues(std::string localInputFilename, float yOffset, std::string normalizationFilename, bool simulatePVFC)
160  {
161  inputFilename = localInputFilename;
162  try
163  {
164  if (verbose) std::cout << "Loading input image..." << std::endl;
165  m_inputImage = mitk::IOUtil::Load<mitk::Image>(inputFilename);
166  if (verbose) std::cout << "Loading input image... [Done]" << std::endl;
167  }
168  catch (...)
169  {
170  if (verbose) std::cout << "No .nrrd file found ... switching to legacy mode." << std::endl;
171  }
172 
173  try
174  {
175  if (simulatePVFC && !normalizationFilename.empty())
176  m_normalizationImage = mitk::IOUtil::Load<mitk::Image>(normalizationFilename);
177  }
178  catch (...)
179  {
180  if (verbose) std::cout << "No normalization .nrrd file found ... will not normalize PVFC" << std::endl;
181  }
182 
183  if (m_normalizationImage.IsNotNull())
184  {
185  mitk::ImageReadAccessor readAccess3(m_normalizationImage, m_normalizationImage->GetVolumeData(0));
186  normalizationVector = (double *)readAccess3.GetData();
187  }
188 
189  if (m_inputImage.IsNotNull()) // load stuff from nrrd file
190  {
191  simulationTimeFromFile = 0;
192 
193  Nx = m_inputImage->GetDimensions()[1];
194  Ny = m_inputImage->GetDimensions()[0];
195  Nz = m_inputImage->GetDimensions()[2];
196 
197  xSpacing = m_inputImage->GetGeometry(0)->GetSpacing()[0];
198  ySpacing = m_inputImage->GetGeometry(0)->GetSpacing()[1];
199  zSpacing = m_inputImage->GetGeometry(0)->GetSpacing()[2];
200 
201  mcflag = std::stoi(m_inputImage->GetProperty("mcflag")->GetValueAsString().c_str()); // mcflag, 0 = uniform, 1 = Gaussian, 2 = iso-pt, 4 = monospectral fraunhofer setup
202  launchflag = std::stoi(m_inputImage->GetProperty("launchflag")->GetValueAsString().c_str());// 0 = let mcxyz calculate trajectory, 1 = manually set launch vector
203  boundaryflag = std::stoi(m_inputImage->GetProperty("boundaryflag")->GetValueAsString().c_str());// 0 = no boundaries, 1 = escape at boundaries, 2 = escape at surface only
204 
205  xs = std::stod(m_inputImage->GetProperty("launchPointX")->GetValueAsString().c_str());
206  ys = std::stod(m_inputImage->GetProperty("launchPointY")->GetValueAsString().c_str()) + yOffset;
207  zs = std::stod(m_inputImage->GetProperty("launchPointZ")->GetValueAsString().c_str());
208 
209  xfocus = std::stod(m_inputImage->GetProperty("focusPointX")->GetValueAsString().c_str());
210  yfocus = std::stod(m_inputImage->GetProperty("focusPointY")->GetValueAsString().c_str());
211  zfocus = std::stod(m_inputImage->GetProperty("focusPointZ")->GetValueAsString().c_str());
212 
213  ux0 = std::stod(m_inputImage->GetProperty("trajectoryVectorX")->GetValueAsString().c_str());
214  uy0 = std::stod(m_inputImage->GetProperty("trajectoryVectorY")->GetValueAsString().c_str());
215  uz0 = std::stod(m_inputImage->GetProperty("trajectoryVectorZ")->GetValueAsString().c_str());
216 
217  radius = std::stod(m_inputImage->GetProperty("radius")->GetValueAsString().c_str());
218  waist = std::stod(m_inputImage->GetProperty("waist")->GetValueAsString().c_str());
219 
220  totalNumberOfVoxels = Nx*Ny*Nz;
221  if (verbose) std::cout << totalNumberOfVoxels << " = sizeof totalNumberOfVoxels" << std::endl;
222 
223  muaVector = (double *)malloc(totalNumberOfVoxels * sizeof(double)); /* tissue structure */
224  musVector = (double *)malloc(totalNumberOfVoxels * sizeof(double)); /* tissue structure */
225  gVector = (double *)malloc(totalNumberOfVoxels * sizeof(double)); /* tissue structure */
226 
227  mitk::ImageReadAccessor readAccess0(m_inputImage, m_inputImage->GetVolumeData(0));
228  muaVector = (double *)readAccess0.GetData();
229  mitk::ImageReadAccessor readAccess1(m_inputImage, m_inputImage->GetVolumeData(1));
230  musVector = (double *)readAccess1.GetData();
231  mitk::ImageReadAccessor readAccess2(m_inputImage, m_inputImage->GetVolumeData(2));
232  gVector = (double *)readAccess2.GetData();
233  }
234  else
235  {
236  mitkThrow() << "No longer support loading of binary tissue files.";
237  }
238  }
239 };
240 
241 class ReturnValues
242 {
243 private:
244  long i1 = 0, i2 = 31; // used Random Generator
245  long ma[56]; // used Random Generator /* ma[0] is not used. */
246  long mj, mk;
247  short i, ii;
248 public:
249  long long Nphotons;
250  double* totalFluence;
251  std::string myname;
252  DetectorVoxel* detectorVoxel;
253 
254  ReturnValues()
255  {
256  detectorVoxel = nullptr;
257  Nphotons = 0;
258  totalFluence = nullptr;
259  }
260 
261  /* SUBROUTINES */
262 
263  /**************************************************************************
264  * RandomGen
265  * A random number generator that generates uniformly
266  * distributed random numbers between 0 and 1 inclusive.
267  * The algorithm is based on:
268  * W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P.
269  * Flannery, "Numerical Recipes in C," Cambridge University
270  * Press, 2nd edition, (1992).
271  * and
272  * D.E. Knuth, "Seminumerical Algorithms," 2nd edition, vol. 2
273  * of "The Art of Computer Programming", Addison-Wesley, (1981).
274  *
275  * When Type is 0, sets Seed as the seed. Make sure 0<Seed<32000.
276  * When Type is 1, returns a random number.
277  * When Type is 2, gets the status of the generator.
278  * When Type is 3, restores the status of the generator.
279  *
280  * The status of the generator is represented by Status[0..56].
281  *
282  * Make sure you initialize the seed before you get random
283  * numbers.
284  ****/
285 #define MBIG 1000000000
286 #define MSEED 161803398
287 #define MZ 0
288 #define FAC 1.0E-9
289 
290  double RandomGen(char Type, long Seed, long *Status) {
291  if (Type == 0) { /* set seed. */
292  if (verbose) std::cout << "Initialized random generator " << this << " with seed: " << Seed << std::endl;
293  mj = MSEED - (Seed < 0 ? -Seed : Seed);
294  mj %= MBIG;
295  ma[55] = mj;
296  mk = 1;
297  for (i = 1; i <= 54; i++) {
298  ii = (21 * i) % 55;
299  ma[ii] = mk;
300  mk = mj - mk;
301  if (mk < MZ)
302  mk += MBIG;
303  mj = ma[ii];
304  }
305  for (ii = 1; ii <= 4; ii++)
306  for (i = 1; i <= 55; i++) {
307  ma[i] -= ma[1 + (i + 30) % 55];
308  if (ma[i] < MZ)
309  ma[i] += MBIG;
310  }
311  i1 = 0;
312  i2 = 31;
313  }
314  else if (Type == 1) { /* get a number. */
315  if (++i1 == 56)
316  i1 = 1;
317  if (++i2 == 56)
318  i2 = 1;
319  mj = ma[i1] - ma[i2];
320  if (mj < MZ)
321  mj += MBIG;
322  ma[i1] = mj;
323  return (mj * FAC);
324  }
325  else if (Type == 2) { /* get status. */
326  for (i = 0; i < 55; i++)
327  Status[i] = ma[i + 1];
328  Status[55] = i1;
329  Status[56] = i2;
330  }
331  else if (Type == 3) { /* restore status. */
332  for (i = 0; i < 55; i++)
333  ma[i + 1] = Status[i];
334  i1 = Status[55];
335  i2 = Status[56];
336  }
337  else
338  puts("Wrong parameter to RandomGen().");
339  return (0);
340  }
341 #undef MBIG
342 #undef MSEED
343 #undef MZ
344 #undef FAC
345 
346  /***********************************************************
347  * Determine if the two position are located in the same voxel
348  * Returns 1 if same voxel, 0 if not same voxel.
349  ****/
350  bool SameVoxel(double x1, double y1, double z1, double x2, double y2, double z2, double dx, double dy, double dz)
351  {
352  double xmin = min2((floor)(x1 / dx), (floor)(x2 / dx))*dx;
353  double ymin = min2((floor)(y1 / dy), (floor)(y2 / dy))*dy;
354  double zmin = min2((floor)(z1 / dz), (floor)(z2 / dz))*dz;
355  double xmax = xmin + dx;
356  double ymax = ymin + dy;
357  double zmax = zmin + dz;
358  bool sv = 0;
359 
360  sv = (x1 <= xmax && x2 <= xmax && y1 <= ymax && y2 <= ymax && z1 < zmax && z2 <= zmax);
361  return (sv);
362  }
363 
364  /***********************************************************
365  * max2
366  ****/
367  double max2(double a, double b) {
368  double m;
369  if (a > b)
370  m = a;
371  else
372  m = b;
373  return m;
374  }
375 
376  /***********************************************************
377  * min2
378  ****/
379  double min2(double a, double b) {
380  double m;
381  if (a >= b)
382  m = b;
383  else
384  m = a;
385  return m;
386  }
387  /***********************************************************
388  * min3
389  ****/
390  double min3(double a, double b, double c) {
391  double m;
392  if (a <= min2(b, c))
393  m = a;
394  else if (b <= min2(a, c))
395  m = b;
396  else
397  m = c;
398  return m;
399  }
400 
401  /********************
402  * my version of FindVoxelFace for no scattering.
403  * s = ls + FindVoxelFace2(x,y,z, tempx, tempy, tempz, dx, dy, dz, ux, uy, uz);
404  ****/
405  double FindVoxelFace2(double x1, double y1, double z1, double /*x2*/, double /*y2*/, double /*z2*/, double dx, double dy, double dz, double ux, double uy, double uz)
406  {
407  int ix1 = floor(x1 / dx);
408  int iy1 = floor(y1 / dy);
409  int iz1 = floor(z1 / dz);
410 
411  int ix2, iy2, iz2;
412  if (ux >= 0)
413  ix2 = ix1 + 1;
414  else
415  ix2 = ix1;
416 
417  if (uy >= 0)
418  iy2 = iy1 + 1;
419  else
420  iy2 = iy1;
421 
422  if (uz >= 0)
423  iz2 = iz1 + 1;
424  else
425  iz2 = iz1;
426 
427  double xs = fabs((ix2*dx - x1) / ux);
428  double ys = fabs((iy2*dy - y1) / uy);
429  double zs = fabs((iz2*dz - z1) / uz);
430 
431  double s = min3(xs, ys, zs);
432 
433  return s;
434  }
435 
436  /***********************************************************
437  * FRESNEL REFLECTANCE
438  * Computes reflectance as photon passes from medium 1 to
439  * medium 2 with refractive indices n1,n2. Incident
440  * angle a1 is specified by cosine value ca1 = cos(a1).
441  * Program returns value of transmitted angle a1 as
442  * value in *ca2_Ptr = cos(a2).
443  ****/
444  double RFresnel(double n1, /* incident refractive index.*/
445  double n2, /* transmit refractive index.*/
446  double ca1, /* cosine of the incident */
447  /* angle a1, 0<a1<90 degrees. */
448  double *ca2_Ptr) /* pointer to the cosine */
449 /* of the transmission */
450 /* angle a2, a2>0. */
451  {
452  double r;
453 
454  if (n1 == n2) {
455  *ca2_Ptr = ca1;
456  r = 0.0;
457  }
458  else if (ca1 > (1.0 - 1.0e-12)) {
459  *ca2_Ptr = ca1;
460  r = (n2 - n1) / (n2 + n1);
461  r *= r;
462  }
463  else if (ca1 < 1.0e-6) {
464  *ca2_Ptr = 0.0;
465  r = 1.0;
466  }
467  else {
468  double sa1, sa2; /* sine of incident and transmission angles. */
469  double ca2; /* cosine of transmission angle. */
470  sa1 = sqrt(1 - ca1*ca1);
471  sa2 = n1*sa1 / n2;
472  if (sa2 >= 1.0) {
473  /* double check for total internal reflection. */
474  *ca2_Ptr = 0.0;
475  r = 1.0;
476  }
477  else {
478  double cap, cam; /* cosines of sum ap or diff am of the two */
479  /* angles: ap = a1 + a2, am = a1 - a2. */
480  double sap, sam; /* sines. */
481  *ca2_Ptr = ca2 = sqrt(1 - sa2*sa2);
482  cap = ca1*ca2 - sa1*sa2; /* c+ = cc - ss. */
483  cam = ca1*ca2 + sa1*sa2; /* c- = cc + ss. */
484  sap = sa1*ca2 + ca1*sa2; /* s+ = sc + cs. */
485  sam = sa1*ca2 - ca1*sa2; /* s- = sc - cs. */
486  r = 0.5*sam*sam*(cam*cam + cap*cap) / (sap*sap*cam*cam);
487  /* rearranged for speed. */
488  }
489  }
490  return(r);
491  } /******** END SUBROUTINE **********/
492 
493  /***********************************************************
494  * the boundary is the face of some voxel
495  * find the the photon's hitting position on the nearest face of the voxel and update the step size.
496  ****/
497  double FindVoxelFace(double x1, double y1, double z1, double x2, double y2, double z2, double dx, double dy, double dz, double ux, double uy, double uz)
498  {
499  double x_1 = x1 / dx;
500  double y_1 = y1 / dy;
501  double z_1 = z1 / dz;
502  double x_2 = x2 / dx;
503  double y_2 = y2 / dy;
504  double z_2 = z2 / dz;
505  double fx_1 = floor(x_1);
506  double fy_1 = floor(y_1);
507  double fz_1 = floor(z_1);
508  double fx_2 = floor(x_2);
509  double fy_2 = floor(y_2);
510  double fz_2 = floor(z_2);
511  double x = 0, y = 0, z = 0, x0 = 0, y0 = 0, z0 = 0, s = 0;
512 
513  if ((fx_1 != fx_2) && (fy_1 != fy_2) && (fz_1 != fz_2)) { //#10
514  fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added
515  fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added
516  fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added
517 
518  x = (max2(fx_1, fx_2) - x_1) / ux;
519  y = (max2(fy_1, fy_2) - y_1) / uy;
520  z = (max2(fz_1, fz_2) - z_1) / uz;
521  if (x == min3(x, y, z)) {
522  x0 = max2(fx_1, fx_2);
523  y0 = (x0 - x_1) / ux*uy + y_1;
524  z0 = (x0 - x_1) / ux*uz + z_1;
525  }
526  else if (y == min3(x, y, z)) {
527  y0 = max2(fy_1, fy_2);
528  x0 = (y0 - y_1) / uy*ux + x_1;
529  z0 = (y0 - y_1) / uy*uz + z_1;
530  }
531  else {
532  z0 = max2(fz_1, fz_2);
533  y0 = (z0 - z_1) / uz*uy + y_1;
534  x0 = (z0 - z_1) / uz*ux + x_1;
535  }
536  }
537  else if ((fx_1 != fx_2) && (fy_1 != fy_2)) { //#2
538  fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added
539  fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added
540  x = (max2(fx_1, fx_2) - x_1) / ux;
541  y = (max2(fy_1, fy_2) - y_1) / uy;
542  if (x == min2(x, y)) {
543  x0 = max2(fx_1, fx_2);
544  y0 = (x0 - x_1) / ux*uy + y_1;
545  z0 = (x0 - x_1) / ux*uz + z_1;
546  }
547  else {
548  y0 = max2(fy_1, fy_2);
549  x0 = (y0 - y_1) / uy*ux + x_1;
550  z0 = (y0 - y_1) / uy*uz + z_1;
551  }
552  }
553  else if ((fy_1 != fy_2) && (fz_1 != fz_2)) { //#3
554  fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added
555  fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added
556  y = (max2(fy_1, fy_2) - y_1) / uy;
557  z = (max2(fz_1, fz_2) - z_1) / uz;
558  if (y == min2(y, z)) {
559  y0 = max2(fy_1, fy_2);
560  x0 = (y0 - y_1) / uy*ux + x_1;
561  z0 = (y0 - y_1) / uy*uz + z_1;
562  }
563  else {
564  z0 = max2(fz_1, fz_2);
565  x0 = (z0 - z_1) / uz*ux + x_1;
566  y0 = (z0 - z_1) / uz*uy + y_1;
567  }
568  }
569  else if ((fx_1 != fx_2) && (fz_1 != fz_2)) { //#4
570  fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added
571  fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added
572  x = (max2(fx_1, fx_2) - x_1) / ux;
573  z = (max2(fz_1, fz_2) - z_1) / uz;
574  if (x == min2(x, z)) {
575  x0 = max2(fx_1, fx_2);
576  y0 = (x0 - x_1) / ux*uy + y_1;
577  z0 = (x0 - x_1) / ux*uz + z_1;
578  }
579  else {
580  z0 = max2(fz_1, fz_2);
581  x0 = (z0 - z_1) / uz*ux + x_1;
582  y0 = (z0 - z_1) / uz*uy + y_1;
583  }
584  }
585  else if (fx_1 != fx_2) { //#5
586  fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added
587  x0 = max2(fx_1, fx_2);
588  y0 = (x0 - x_1) / ux*uy + y_1;
589  z0 = (x0 - x_1) / ux*uz + z_1;
590  }
591  else if (fy_1 != fy_2) { //#6
592  fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added
593  y0 = max2(fy_1, fy_2);
594  x0 = (y0 - y_1) / uy*ux + x_1;
595  z0 = (y0 - y_1) / uy*uz + z_1;
596  }
597  else { //#7
598  z0 = max2(fz_1, fz_2);
599  fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added
600  x0 = (z0 - z_1) / uz*ux + x_1;
601  y0 = (z0 - z_1) / uz*uy + y_1;
602  }
603  //s = ( (x0-fx_1)*dx + (y0-fy_1)*dy + (z0-fz_1)*dz )/3;
604  //s = sqrt( SQR((x0-x_1)*dx) + SQR((y0-y_1)*dy) + SQR((z0-z_1)*dz) );
605  //s = sqrt(SQR(x0-x_1)*SQR(dx) + SQR(y0-y_1)*SQR(dy) + SQR(z0-z_1)*SQR(dz));
606  s = sqrt(SQR((x0 - x_1)*dx) + SQR((y0 - y_1)*dy) + SQR((z0 - z_1)*dz));
607  return (s);
608  }
609 };
610 
611 /* DECLARE FUNCTIONS */
612 
613 void runMonteCarlo(InputValues* inputValues, ReturnValues* returnValue, int thread, mitk::pa::MonteCarloThreadHandler::Pointer threadHandler);
614 
615 int detector_x = -1;
616 int detector_z = -1;
617 bool interpretAsTime = true;
618 bool simulatePVFC = false;
620 float requestedSimulationTime = 0; // in minutes
622 float yOffset = 0; // in mm
623 bool saveLegacy = false;
625 std::string inputFilename;
626 std::string outputFilename;
627 
628 mitk::pa::Probe::Pointer m_PhotoacousticProbe;
629 
630 int main(int argc, char * argv[]) {
631  mitkCommandLineParser parser;
632  // set general information
633  parser.setCategory("MITK-Photoacoustics");
634  parser.setTitle("Mitk MCxyz");
635  parser.setDescription("Runs Monte Carlo simulations on inputed tissues.");
636  parser.setContributor("CAI, DKFZ based on code by Jacques and Li");
637 
638  // how should arguments be prefixed
639  parser.setArgumentPrefix("--", "-");
640  // add each argument, unless specified otherwise each argument is optional
641  // see mitkCommandLineParser::addArgument for more information
642  parser.beginGroup("Required I/O parameters");
643  parser.addArgument(
644  "input", "i", mitkCommandLineParser::File,
645  "Input tissue file", "input tissue file (*.nrrd)",
646  us::Any(), false, false, false, mitkCommandLineParser::Input);
647  parser.addArgument(
648  "output", "o", mitkCommandLineParser::File,
649  "Output fluence file", "where to save the simulated fluence (*.nrrd)",
650  us::Any(), false, false, false, mitkCommandLineParser::Output);
651  parser.endGroup();
652  parser.beginGroup("Optional parameters");
653  parser.addArgument(
654  "verbose", "v", mitkCommandLineParser::Bool,
655  "Verbose Output", "Whether to produce verbose, or rather debug output");
656  parser.addArgument(
657  "detector-x", "dx", mitkCommandLineParser::Int,
658  "Detector voxel x position", "Determines the x position of the detector voxel (default: -1 = dont use detector voxel)", -1);
659  parser.addArgument(
660  "detector-z", "dz", mitkCommandLineParser::Int,
661  "Detector voxel z position", "Determines the z position of the detector voxel (default: -1 = dont use detector voxel)", -1);
662  parser.addArgument(
663  "number-of-photons", "n", mitkCommandLineParser::Int,
664  "Number of photons", "Specifies the number of photons (default: 100000). Simulation stops after that number. Use -t --timer to define a timer instead");
665  parser.addArgument(
666  "timer", "t", mitkCommandLineParser::Float,
667  "Simulation time in min", "Specifies the amount of time for simutation (default: 0). Simulation stops after that number of minutes. -n --number-of-photons is the override and default behavior and defines the maximum number of photons instead. If no simulation time or number of photons is specified the file time is taken.");
668  parser.addArgument(
669  "y-offset", "yo", mitkCommandLineParser::Float,
670  "Probe Y-Offset in mm", "Specifies an offset of the photoacoustic probe in the y direction depending on the initial probe position (default: 0) in mm.");
671  parser.addArgument(
672  "jobs", "j", mitkCommandLineParser::Int,
673  "Number of jobs", "Specifies the number of jobs for simutation (default: -1 which starts as many jobs as supported).");
674  parser.addArgument(
675  "probe-xml", "p", mitkCommandLineParser::File,
676  "Xml definition of the probe", "Specifies the absolute path of the location of the xml definition file of the probe design.", us::Any(), true, false, false, mitkCommandLineParser::Input);
677  parser.addArgument("normalization-file", "nf", mitkCommandLineParser::File,
678  "Input normalization file", "The input normalization file is used for normalization of the number of photons in the PVFC calculations.", us::Any(), true, false, false, mitkCommandLineParser::Input);
679  parser.endGroup();
680 
681  // parse arguments, this method returns a mapping of long argument names and their values
682  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
683  if (parsedArgs.size() == 0)
684  return EXIT_FAILURE;
685  // parse, cast and set required arguments
686  inputFilename = us::any_cast<std::string>(parsedArgs["input"]);
687  // strip ending
688  inputFilename = inputFilename.substr(0, inputFilename.find("_H.mci"));
689  inputFilename = inputFilename.substr(0, inputFilename.find("_T.bin"));
690 
691  outputFilename = us::any_cast<std::string>(parsedArgs["output"]);
692  // add .nrrd if not there
693  std::string suffix = ".nrrd";
694  if (outputFilename.compare(outputFilename.size() - suffix.size(), suffix.size(), suffix) != 0)
695  outputFilename = outputFilename + suffix;
696 
697  // default values for optional arguments
698  // parse, cast and set optional arguments if given
699  if (parsedArgs.count("verbose"))
700  {
701  verbose = us::any_cast<bool>(parsedArgs["verbose"]);
702  }
703  if (parsedArgs.count("detector-x"))
704  {
705  detector_x = us::any_cast<int>(parsedArgs["detector-x"]);
706  }
707  if (parsedArgs.count("detector-z"))
708  {
709  detector_z = us::any_cast<int>(parsedArgs["detector-z"]);
710  }
711  if (parsedArgs.count("timer"))
712  {
713  requestedSimulationTime = us::any_cast<float>(parsedArgs["timer"]);
715  }
716  if (parsedArgs.count("y-offset"))
717  {
718  yOffset = us::any_cast<float>(parsedArgs["y-offset"]);
719  }
720  if (parsedArgs.count("number-of-photons"))
721  {
722  requestedNumberOfPhotons = us::any_cast<int>(parsedArgs["number-of-photons"]);
723  if (requestedNumberOfPhotons > 0) interpretAsTime = false;
724  }
725  if (parsedArgs.count("jobs"))
726  {
727  concurentThreadsSupported = us::any_cast<int>(parsedArgs["jobs"]);
728  }
729  if (parsedArgs.count("probe-xml"))
730  {
731  std::string inputXmlProbeDesign = us::any_cast<std::string>(parsedArgs["probe-xml"]);
732  m_PhotoacousticProbe = mitk::pa::Probe::New(inputXmlProbeDesign, verbose);
733  if (!m_PhotoacousticProbe->IsValid())
734  {
735  std::cerr << "Xml File was not valid. Simulation failed." << std::endl;
736  return EXIT_FAILURE;
737  }
738  }
739  if (parsedArgs.count("normalization-file"))
740  {
741  normalizationFilename = us::any_cast<std::string>(parsedArgs["normalization-file"]);
742  }
743 
745  {
746  concurentThreadsSupported = std::thread::hardware_concurrency();
747  if (concurentThreadsSupported == 0)
748  {
749  std::cout << "Could not determine number of available threads. Launching only one." << std::endl;
751  }
752  }
753 
754  if (detector_x != -1 && detector_z != -1)
755  {
756  if (verbose)
757  std::cout << "Performing PVFC calculation for x=" << detector_x << " and z=" << detector_z << std::endl;
758  simulatePVFC = true;
759  }
760  else
761  {
762  if (verbose)
763  std::cout << "Will not perform PVFC calculation due to x=" << detector_x << " and/or z=" << detector_z << std::endl;
764  }
765 
766  InputValues allInput = InputValues();
767  allInput.LoadValues(inputFilename, yOffset, normalizationFilename, simulatePVFC);
768 
769  std::vector<ReturnValues> allValues(concurentThreadsSupported);
770  auto* threads = new std::thread[concurentThreadsSupported];
771 
772  for (long i = 0; i < concurentThreadsSupported; i++)
773  {
774  auto* tmp = new ReturnValues();
775  allValues.push_back(*tmp);
776  }
777 
778  if (verbose) std::cout << "Initializing MonteCarloThreadHandler" << std::endl;
779 
780  long timeMetric;
781  if (interpretAsTime)
782  {
784  requestedSimulationTime = allInput.simulationTimeFromFile;
785 
786  timeMetric = requestedSimulationTime * 60 * 1000;
787  }
788  else
789  {
790  timeMetric = requestedNumberOfPhotons;
791  }
792 
793  mitk::pa::MonteCarloThreadHandler::Pointer threadHandler = mitk::pa::MonteCarloThreadHandler::New(timeMetric, interpretAsTime);
794 
795  if (simulatePVFC)
796  threadHandler->SetPackageSize(1000);
797 
798  if (verbose) std::cout << "\nStarting simulation ...\n" << std::endl;
799 
800  auto simulationStartTime = std::chrono::system_clock::now();
801 
802  for (int i = 0; i < concurentThreadsSupported; i++)
803  {
804  threads[i] = std::thread(runMonteCarlo, &allInput, &allValues[i], (i + 1), threadHandler);
805  }
806 
807  for (int i = 0; i < concurentThreadsSupported; i++)
808  {
809  threads[i].join();
810  }
811 
812  auto simulationFinishTime = std::chrono::system_clock::now();
813  auto simulationTimeElapsed = simulationFinishTime - simulationStartTime;
814 
815  if (verbose) std::cout << "\n\nFinished simulation\n\n" << std::endl;
816  std::cout << "total time for simulation: "
817  << (int)std::chrono::duration_cast<std::chrono::seconds>(simulationTimeElapsed).count() << "sec " << std::endl;
818 
819  /**** SAVE
820  Convert data to relative fluence rate [cm^-2] and save.
821  *****/
822 
823  if (!simulatePVFC)
824  {
825  if (verbose) std::cout << "Allocating memory for normal simulation result ... ";
826  auto* finalTotalFluence = (double *)malloc(allInput.totalNumberOfVoxels * sizeof(double));
827  if (verbose) std::cout << "[OK]" << std::endl;
828  if (verbose) std::cout << "Cleaning memory for normal simulation result ...";
829  for (int i = 0; i < allInput.totalNumberOfVoxels; i++) {
830  finalTotalFluence[i] = 0;
831  }
832  if (verbose) std::cout << "[OK]" << std::endl;
833 
834  if (verbose) std::cout << "Calculating resulting fluence ... ";
835  double tdx = 0, tdy = 0, tdz = 0;
836  long long tNphotons = 0;
837  for (int t = 0; t < concurentThreadsSupported; t++)
838  {
839  tdx = allInput.xSpacing;
840  tdy = allInput.ySpacing;
841  tdz = allInput.zSpacing;
842  tNphotons += allValues[t].Nphotons;
843  for (int voxelNumber = 0; voxelNumber < allInput.totalNumberOfVoxels; voxelNumber++) {
844  finalTotalFluence[voxelNumber] += allValues[t].totalFluence[voxelNumber];
845  }
846  }
847  if (verbose) std::cout << "[OK]" << std::endl;
848  std::cout << "total number of photons simulated: "
849  << tNphotons << std::endl;
850 
851  // Normalize deposition (A) to yield fluence rate (F).
852  double temp = tdx*tdy*tdz*tNphotons;
853  for (int i = 0; i < allInput.totalNumberOfVoxels; i++) {
854  finalTotalFluence[i] /= temp*allInput.muaVector[i];
855  }
856 
857  if (verbose) std::cout << "Saving normal simulated fluence result to " << outputFilename << " ... ";
858 
859  mitk::Image::Pointer resultImage = mitk::Image::New();
860  mitk::PixelType TPixel = mitk::MakeScalarPixelType<double>();
861  auto* dimensionsOfImage = new unsigned int[3];
862 
863  // Copy dimensions
864  dimensionsOfImage[0] = allInput.Ny;
865  dimensionsOfImage[1] = allInput.Nx;
866  dimensionsOfImage[2] = allInput.Nz;
867 
868  resultImage->Initialize(TPixel, 3, dimensionsOfImage);
869 
870  mitk::Vector3D spacing;
871  spacing[0] = allInput.ySpacing;
872  spacing[1] = allInput.xSpacing;
873  spacing[2] = allInput.zSpacing;
874  resultImage->SetSpacing(spacing);
875  resultImage->SetImportVolume(finalTotalFluence, 0, 0, mitk::Image::CopyMemory);
876 
877  resultImage->GetPropertyList()->SetFloatProperty("y-offset", yOffset);
879 
880  mitk::IOUtil::Save(resultImage, outputFilename);
881 
882  if (verbose) std::cout << "[OK]" << std::endl;
883 
884  if (verbose)
885  {
886  std::cout << "x spacing = " << tdx << std::endl;
887  std::cout << "y spacing = " << tdy << std::endl;
888  std::cout << "z spacing = " << tdz << std::endl;
889  std::cout << "total number of voxels = " << allInput.totalNumberOfVoxels << std::endl;
890  std::cout << "number of photons = " << (int)tNphotons << std::endl;
891  }
892  }
893  else // if simulate PVFC
894  {
895  if (verbose) std::cout << "Allocating memory for PVFC simulation result ... ";
896  double* detectorFluence = ((double*)malloc(allInput.totalNumberOfVoxels * sizeof(double)));
897  if (verbose) std::cout << "[OK]" << std::endl;
898  if (verbose) std::cout << "Cleaning memory for PVFC simulation result ...";
899  for (int i = 0; i < allInput.totalNumberOfVoxels; i++) {
900  detectorFluence[i] = 0;
901  }
902  if (verbose) std::cout << "[OK]" << std::endl;
903 
904  if (verbose) std::cout << "Calculating resulting PVFC fluence ... ";
905  double tdx = 0, tdy = 0, tdz = 0;
906  long long tNphotons = 0;
907  long pvfcPhotons = 0;
908  for (int t = 0; t < concurentThreadsSupported; t++)
909  {
910  tdx = allInput.xSpacing;
911  tdy = allInput.ySpacing;
912  tdz = allInput.zSpacing;
913  tNphotons += allValues[t].Nphotons;
914  pvfcPhotons += allValues[t].detectorVoxel->m_NumberPhotonsCurrent;
915  for (int voxelNumber = 0; voxelNumber < allInput.totalNumberOfVoxels; voxelNumber++) {
916  detectorFluence[voxelNumber] +=
917  allValues[t].detectorVoxel->fluenceContribution[voxelNumber];
918  }
919  }
920  if (verbose) std::cout << "[OK]" << std::endl;
921  std::cout << "total number of photons simulated: "
922  << tNphotons << std::endl;
923 
924  // Normalize deposition (A) to yield fluence rate (F).
925  double temp = tdx*tdy*tdz*tNphotons;
926  for (int i = 0; i < allInput.totalNumberOfVoxels; i++) {
927  detectorFluence[i] /= temp*allInput.muaVector[i];
928  }
929 
930  if (verbose) std::cout << "Saving PVFC ...";
931 
932  std::stringstream detectorname("");
933  double detectorX = allValues[0].detectorVoxel->location.x;
934  double detectorY = allValues[0].detectorVoxel->location.y;
935  double detectorZ = allValues[0].detectorVoxel->location.z;
936  detectorname << detectorX << "," << detectorY << ","
937  << detectorZ << "FluenceContribution.nrrd";
938  // Save the binary file
939  std::string outputFileBase = outputFilename.substr(0, outputFilename.find(".nrrd"));
940  outputFilename = outputFileBase + "_p" + detectorname.str().c_str();
941 
943  auto* dimensionsOfPvfcImage = new unsigned int[3];
944 
945  // Copy dimensions
946  dimensionsOfPvfcImage[0] = allInput.Ny;
947  dimensionsOfPvfcImage[1] = allInput.Nx;
948  dimensionsOfPvfcImage[2] = allInput.Nz;
949 
950  pvfcImage->Initialize(mitk::MakeScalarPixelType<double>(), 3, dimensionsOfPvfcImage);
951 
952  mitk::Vector3D pvfcSpacing;
953  pvfcSpacing[0] = allInput.ySpacing;
954  pvfcSpacing[1] = allInput.xSpacing;
955  pvfcSpacing[2] = allInput.zSpacing;
956  pvfcImage->SetSpacing(pvfcSpacing);
957  pvfcImage->SetImportVolume(detectorFluence, 0, 0, mitk::Image::CopyMemory);
958 
959  pvfcImage->GetPropertyList()->SetFloatProperty("detector-x", detectorX);
961  pvfcImage->GetPropertyList()->SetFloatProperty("detector-y", detectorY);
963  pvfcImage->GetPropertyList()->SetFloatProperty("detector-z", detectorZ);
965  pvfcImage->GetPropertyList()->SetFloatProperty("normalization-factor", allValues[0].detectorVoxel->m_PhotonNormalizationValue);
967  pvfcImage->GetPropertyList()->SetFloatProperty("simulated-photons", pvfcPhotons);
969 
971 
972  if (verbose) std::cout << "[OK]" << std::endl;
973 
974  if (verbose)
975  {
976  std::cout << "x spacing = " << tdx << std::endl;
977  std::cout << "y spacing = " << tdy << std::endl;
978  std::cout << "z spacing = " << tdz << std::endl;
979  std::cout << "total number of voxels = " << allInput.totalNumberOfVoxels << std::endl;
980  std::cout << "number of photons = " << (int)tNphotons << std::endl;
981  }
982  }
983 
984  exit(EXIT_SUCCESS);
985 } /* end of main */
986 
987 /* CORE FUNCTION */
988 void runMonteCarlo(InputValues* inputValues, ReturnValues* returnValue, int thread, mitk::pa::MonteCarloThreadHandler::Pointer threadHandler)
989 {
990  if (verbose) std::cout << "Thread " << thread << ": Locking Mutex ..." << std::endl;
991  if (verbose) std::cout << "[OK]" << std::endl;
992  if (verbose) std::cout << "Initializing ... ";
993 
994  /* Propagation parameters */
995  double x, y, z; /* photon position */
996  double ux, uy, uz; /* photon trajectory as cosines */
997  double uxx, uyy, uzz; /* temporary values used during SPIN */
998  double s; /* step sizes. s = -log(RND)/mus [cm] */
999  double sleft; /* dimensionless */
1000  double costheta; /* cos(theta) */
1001  double sintheta; /* sin(theta) */
1002  double cospsi; /* cos(psi) */
1003  double sinpsi; /* sin(psi) */
1004  double psi; /* azimuthal angle */
1005  long photonIterator = 0; /* current photon */
1006  double W; /* photon weight */
1007  double absorb; /* weighted deposited in a step due to absorption */
1008  short photon_status; /* flag = ALIVE=1 or DEAD=0 */
1009  bool sv; /* Are they in the same voxel? */
1010 
1011  /* dummy variables */
1012  double rnd; /* assigned random value 0-1 */
1013  double r, phi; /* dummy values */
1014  long i, j; /* dummy indices */
1015  double tempx, tempy, tempz; /* temporary variables, used during photon step. */
1016  int ix, iy, iz; /* Added. Used to track photons */
1017  double temp; /* dummy variable */
1018  int bflag; /* boundary flag: 0 = photon inside volume. 1 = outside volume */
1019  int CNT = 0;
1020 
1021  returnValue->totalFluence = (double *)malloc(inputValues->totalNumberOfVoxels * sizeof(double)); /* relative fluence rate [W/cm^2/W.delivered] */
1022 
1023  if (detector_x != -1 && detector_z != -1)
1024  {
1025  if (detector_x<0 || detector_x>inputValues->Nx)
1026  {
1027  std::cout << "Requested detector x position not valid. Needs to be >= 0 and <= " << inputValues->Nx << std::endl;
1028  exit(EXIT_FAILURE);
1029  }
1030  if (detector_z<1 || detector_z>inputValues->Nz)
1031  {
1032  std::cout << "Requested detector z position not valid. Needs to be > 0 and <= " << inputValues->Nz << std::endl;
1033  exit(EXIT_FAILURE);
1034  }
1035 
1036  double photonNormalizationValue = 1 / inputValues->GetNormalizationValue(detector_x, inputValues->Ny / 2, detector_z);
1037  returnValue->detectorVoxel = new DetectorVoxel(initLocation(detector_x, inputValues->Ny / 2, detector_z, 0), inputValues->totalNumberOfVoxels, photonNormalizationValue);
1038  }
1039 
1040  /**** ======================== MAJOR CYCLE ============================ *****/
1041 
1042  auto duration = std::chrono::system_clock::now().time_since_epoch();
1043  returnValue->RandomGen(0, (std::chrono::duration_cast<std::chrono::milliseconds>(duration).count() + thread) % 32000, nullptr); /* initiate with seed = 1, or any long integer. */
1044  for (j = 0; j < inputValues->totalNumberOfVoxels; j++) returnValue->totalFluence[j] = 0; // ensure F[] starts empty.
1045 
1046  /**** RUN Launch N photons, initializing each one before progation. *****/
1047 
1048  long photonsToSimulate = 0;
1049 
1050  do {
1051  photonsToSimulate = threadHandler->GetNextWorkPackage();
1052  if (returnValue->detectorVoxel != nullptr)
1053  {
1054  photonsToSimulate = photonsToSimulate * returnValue->detectorVoxel->m_PhotonNormalizationValue;
1055  }
1056 
1057  if (verbose)
1058  MITK_INFO << "Photons to simulate: " << photonsToSimulate;
1059 
1060  photonIterator = 0L;
1061 
1062  do {
1063  /**** LAUNCH Initialize photon position and trajectory. *****/
1064 
1065  photonIterator += 1; /* increment photon count */
1066  W = 1.0; /* set photon weight to one */
1067  photon_status = ALIVE; /* Launch an ALIVE photon */
1068  CNT = 0;
1069 
1070  /**** SET SOURCE* Launch collimated beam at x,y center.****/
1071  /****************************/
1072  /* Initial position. */
1073 
1074  if (m_PhotoacousticProbe.IsNotNull())
1075  {
1076  double rnd1 = -1;
1077  double rnd2 = -1;
1078  double rnd3 = -1;
1079  double rnd4 = -1;
1080  double rnd5 = -1;
1081  double rnd6 = -1;
1082  double rnd7 = -1;
1083  double rnd8 = -1;
1084 
1085  while ((rnd1 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1086  while ((rnd2 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1087  while ((rnd3 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1088  while ((rnd4 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1089  while ((rnd5 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1090  while ((rnd6 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1091  while ((rnd7 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1092  while ((rnd8 = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1093 
1094  mitk::pa::LightSource::PhotonInformation info = m_PhotoacousticProbe->GetNextPhoton(rnd1, rnd2, rnd3, rnd4, rnd5, rnd6, rnd7, rnd8);
1095  x = info.xPosition;
1096  y = yOffset + info.yPosition;
1097  z = info.zPosition;
1098  ux = info.xAngle;
1099  uy = info.yAngle;
1100  uz = info.zAngle;
1101  if (verbose)
1102  std::cout << "Created photon at position (" << x << "|" << y << "|" << z << ") with angles (" << ux << "|" << uy << "|" << uz << ")." << std::endl;
1103  }
1104  else
1105  {
1106  /* trajectory */
1107  if (inputValues->launchflag == 1) // manually set launch
1108  {
1109  x = inputValues->xs;
1110  y = inputValues->ys;
1111  z = inputValues->zs;
1112  ux = inputValues->ux0;
1113  uy = inputValues->uy0;
1114  uz = inputValues->uz0;
1115  }
1116  else // use mcflag
1117  {
1118  if (inputValues->mcflag == 0) // uniform beam
1119  {
1120  // set launch point and width of beam
1121  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0); // avoids rnd = 0
1122  r = inputValues->radius*sqrt(rnd); // radius of beam at launch point
1123  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0); // avoids rnd = 0
1124  phi = rnd*2.0*PI;
1125  x = inputValues->xs + r*cos(phi);
1126  y = inputValues->ys + r*sin(phi);
1127  z = inputValues->zs;
1128  // set trajectory toward focus
1129  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0); // avoids rnd = 0
1130  r = inputValues->waist*sqrt(rnd); // radius of beam at focus
1131  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0); // avoids rnd = 0
1132  phi = rnd*2.0*PI;
1133 
1134  // !!!!!!!!!!!!!!!!!!!!!!! setting input values will braek
1135 
1136  inputValues->xfocus = r*cos(phi);
1137  inputValues->yfocus = r*sin(phi);
1138  temp = sqrt((x - inputValues->xfocus)*(x - inputValues->xfocus)
1139  + (y - inputValues->yfocus)*(y - inputValues->yfocus) + inputValues->zfocus*inputValues->zfocus);
1140  ux = -(x - inputValues->xfocus) / temp;
1141  uy = -(y - inputValues->yfocus) / temp;
1142  uz = sqrt(1 - ux*ux + uy*uy);
1143  }
1144  else if (inputValues->mcflag == 5) // Multispectral DKFZ prototype
1145  {
1146  // set launch point and width of beam
1147  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1148 
1149  //offset in x direction in cm (random)
1150  x = (rnd*2.5) - 1.25;
1151 
1152  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1153  double b = ((rnd)-0.5);
1154  y = (b > 0 ? yOffset + 1.5 : yOffset - 1.5);
1155  z = 0.1;
1156  ux = 0;
1157 
1158  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1159 
1160  //Angle of beam in y direction
1161  uy = sin((rnd*0.42) - 0.21 + (b < 0 ? 1.0 : -1.0) * 0.436);
1162 
1163  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1164 
1165  // angle of beam in x direction
1166  ux = sin((rnd*0.42) - 0.21);
1167  uz = sqrt(1 - ux*ux - uy*uy);
1168  }
1169  else if (inputValues->mcflag == 4) // Monospectral prototype DKFZ
1170  {
1171  // set launch point and width of beam
1172  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1173 
1174  //offset in x direction in cm (random)
1175  x = (rnd*2.5) - 1.25;
1176 
1177  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1178  double b = ((rnd)-0.5);
1179  y = (b > 0 ? yOffset + 0.83 : yOffset - 0.83);
1180  z = 0.1;
1181  ux = 0;
1182 
1183  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1184 
1185  //Angle of beam in y direction
1186  uy = sin((rnd*0.42) - 0.21 + (b < 0 ? 1.0 : -1.0) * 0.375);
1187 
1188  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1189 
1190  // angle of beam in x direction
1191  ux = sin((rnd*0.42) - 0.21);
1192  uz = sqrt(1 - ux*ux - uy*uy);
1193  }
1194  else { // isotropic pt source
1195  costheta = 1.0 - 2.0 * returnValue->RandomGen(1, 0, nullptr);
1196  sintheta = sqrt(1.0 - costheta*costheta);
1197  psi = 2.0 * PI * returnValue->RandomGen(1, 0, nullptr);
1198  cospsi = cos(psi);
1199  if (psi < PI)
1200  sinpsi = sqrt(1.0 - cospsi*cospsi);
1201  else
1202  sinpsi = -sqrt(1.0 - cospsi*cospsi);
1203  x = inputValues->xs;
1204  y = inputValues->ys;
1205  z = inputValues->zs;
1206  ux = sintheta*cospsi;
1207  uy = sintheta*sinpsi;
1208  uz = costheta;
1209  }
1210  } // end use mcflag
1211  }
1212  /****************************/
1213 
1214  /* Get tissue voxel properties of launchpoint.
1215  * If photon beyond outer edge of defined voxels,
1216  * the tissue equals properties of outermost voxels.
1217  * Therefore, set outermost voxels to infinite background value.
1218  */
1219  ix = (int)(inputValues->Nx / 2 + x / inputValues->xSpacing);
1220  iy = (int)(inputValues->Ny / 2 + y / inputValues->ySpacing);
1221  iz = (int)(z / inputValues->zSpacing);
1222  if (ix >= inputValues->Nx) ix = inputValues->Nx - 1;
1223  if (iy >= inputValues->Ny) iy = inputValues->Ny - 1;
1224  if (iz >= inputValues->Nz) iz = inputValues->Nz - 1;
1225  if (ix < 0) ix = 0;
1226  if (iy < 0) iy = 0;
1227  if (iz < 0) iz = 0;
1228  /* Get the tissue type of located voxel */
1229  i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy);
1230 
1231  bflag = 1; // initialize as 1 = inside volume, but later check as photon propagates.
1232 
1233  if (returnValue->detectorVoxel != nullptr)
1234  returnValue->detectorVoxel->recordedPhotonRoute->clear();
1235 
1236  /* HOP_DROP_SPIN_CHECK
1237  Propagate one photon until it dies as determined by ROULETTE.
1238  *******/
1239  do {
1240  /**** HOP
1241  Take step to new position
1242  s = dimensionless stepsize
1243  x, uy, uz are cosines of current photon trajectory
1244  *****/
1245  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0); /* yields 0 < rnd <= 1 */
1246  sleft = -log(rnd); /* dimensionless step */
1247  CNT += 1;
1248 
1249  do { // while sleft>0
1250  s = sleft / inputValues->musVector[i]; /* Step size [cm].*/
1251  tempx = x + s*ux; /* Update positions. [cm] */
1252  tempy = y + s*uy;
1253  tempz = z + s*uz;
1254 
1255  sv = returnValue->SameVoxel(x, y, z, tempx, tempy, tempz, inputValues->xSpacing, inputValues->ySpacing, inputValues->zSpacing);
1256  if (sv) /* photon in same voxel */
1257  {
1258  x = tempx; /* Update positions. */
1259  y = tempy;
1260  z = tempz;
1261 
1262  /**** DROP
1263  Drop photon weight (W) into local bin.
1264  *****/
1265 
1266  absorb = W*(1 - exp(-inputValues->muaVector[i] * s)); /* photon weight absorbed at this step */
1267  W -= absorb; /* decrement WEIGHT by amount absorbed */
1268  // If photon within volume of heterogeneity, deposit energy in F[].
1269  // Normalize F[] later, when save output.
1270  if (bflag)
1271  {
1272  i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy);
1273  returnValue->totalFluence[i] += absorb;
1274  // only save data if blag==1, i.e., photon inside simulation cube
1275 
1276  //For each detectorvoxel
1277  if (returnValue->detectorVoxel != nullptr)
1278  {
1279  //Add photon position to the recorded photon route
1280  returnValue->detectorVoxel->recordedPhotonRoute->push_back(initLocation(ix, iy, iz, absorb));
1281 
1282  //If the photon is currently at the detector position
1283  if ((returnValue->detectorVoxel->location.x == ix)
1284  && ((returnValue->detectorVoxel->location.y == iy)
1285  || (returnValue->detectorVoxel->location.y - 1 == iy))
1286  && (returnValue->detectorVoxel->location.z == iz))
1287  {
1288  //For each voxel in the recorded photon route
1289  for (unsigned int routeIndex = 0; routeIndex < returnValue->detectorVoxel->recordedPhotonRoute->size(); routeIndex++)
1290  {
1291  //increment the fluence contribution at that particular position
1292  i = (long)(returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).z*inputValues->Ny*inputValues->Nx
1293  + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).x*inputValues->Ny
1294  + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).y);
1295  returnValue->detectorVoxel->fluenceContribution[i] += returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).absorb;
1296  }
1297 
1298  //Clear the recorded photon route
1299  returnValue->detectorVoxel->m_NumberPhotonsCurrent++;
1300  returnValue->detectorVoxel->recordedPhotonRoute->clear();
1301  }
1302  }
1303  }
1304 
1305  /* Update sleft */
1306  sleft = 0; /* dimensionless step remaining */
1307  }
1308  else /* photon has crossed voxel boundary */
1309  {
1310  /* step to voxel face + "littlest step" so just inside new voxel. */
1311  s = ls + returnValue->FindVoxelFace2(x, y, z, tempx, tempy, tempz, inputValues->xSpacing, inputValues->ySpacing, inputValues->zSpacing, ux, uy, uz);
1312 
1313  /**** DROP
1314  Drop photon weight (W) into local bin.
1315  *****/
1316  absorb = W*(1 - exp(-inputValues->muaVector[i] * s)); /* photon weight absorbed at this step */
1317  W -= absorb; /* decrement WEIGHT by amount absorbed */
1318  // If photon within volume of heterogeneity, deposit energy in F[].
1319  // Normalize F[] later, when save output.
1320  if (bflag)
1321  {
1322  // only save data if bflag==1, i.e., photon inside simulation cube
1323 
1324  //For each detectorvoxel
1325  if (returnValue->detectorVoxel != nullptr)
1326  {
1327  //Add photon position to the recorded photon route
1328  returnValue->detectorVoxel->recordedPhotonRoute->push_back(initLocation(ix, iy, iz, absorb));
1329 
1330  //If the photon is currently at the detector position
1331  if ((returnValue->detectorVoxel->location.x == ix)
1332  && ((returnValue->detectorVoxel->location.y == iy)
1333  || (returnValue->detectorVoxel->location.y - 1 == iy))
1334  && (returnValue->detectorVoxel->location.z == iz))
1335  {
1336  //For each voxel in the recorded photon route
1337  for (unsigned int routeIndex = 0; routeIndex < returnValue->detectorVoxel->recordedPhotonRoute->size(); routeIndex++)
1338  {
1339  //increment the fluence contribution at that particular position
1340  i = (long)(returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).z*inputValues->Ny*inputValues->Nx
1341  + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).x*inputValues->Ny
1342  + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).y);
1343  returnValue->detectorVoxel->fluenceContribution[i] += returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).absorb;
1344  }
1345 
1346  //Clear the recorded photon route
1347  returnValue->detectorVoxel->m_NumberPhotonsCurrent++;
1348  returnValue->detectorVoxel->recordedPhotonRoute->clear();
1349  }
1350  }
1351 
1352  i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy);
1353  returnValue->totalFluence[i] += absorb;
1354  }
1355 
1356  /* Update sleft */
1357  sleft -= s*inputValues->musVector[i]; /* dimensionless step remaining */
1358  if (sleft <= ls) sleft = 0;
1359 
1360  /* Update positions. */
1361  x += s*ux;
1362  y += s*uy;
1363  z += s*uz;
1364 
1365  // pointers to voxel containing optical properties
1366  ix = (int)(inputValues->Nx / 2 + x / inputValues->xSpacing);
1367  iy = (int)(inputValues->Ny / 2 + y / inputValues->ySpacing);
1368  iz = (int)(z / inputValues->zSpacing);
1369 
1370  bflag = 1; // Boundary flag. Initialize as 1 = inside volume, then check.
1371  if (inputValues->boundaryflag == 0) { // Infinite medium.
1372  // Check if photon has wandered outside volume.
1373  // If so, set tissue type to boundary value, but let photon wander.
1374  // Set blag to zero, so DROP does not deposit energy.
1375  if (iz >= inputValues->Nz) { iz = inputValues->Nz - 1; bflag = 0; }
1376  if (ix >= inputValues->Nx) { ix = inputValues->Nx - 1; bflag = 0; }
1377  if (iy >= inputValues->Ny) { iy = inputValues->Ny - 1; bflag = 0; }
1378  if (iz < 0) { iz = 0; bflag = 0; }
1379  if (ix < 0) { ix = 0; bflag = 0; }
1380  if (iy < 0) { iy = 0; bflag = 0; }
1381  }
1382  else if (inputValues->boundaryflag == 1) { // Escape at boundaries
1383  if (iz >= inputValues->Nz) { iz = inputValues->Nz - 1; photon_status = DEAD; sleft = 0; }
1384  if (ix >= inputValues->Nx) { ix = inputValues->Nx - 1; photon_status = DEAD; sleft = 0; }
1385  if (iy >= inputValues->Ny) { iy = inputValues->Ny - 1; photon_status = DEAD; sleft = 0; }
1386  if (iz < 0) { iz = 0; photon_status = DEAD; sleft = 0; }
1387  if (ix < 0) { ix = 0; photon_status = DEAD; sleft = 0; }
1388  if (iy < 0) { iy = 0; photon_status = DEAD; sleft = 0; }
1389  }
1390  else if (inputValues->boundaryflag == 2) { // Escape at top surface, no x,y bottom z boundaries
1391  if (iz >= inputValues->Nz) { iz = inputValues->Nz - 1; bflag = 0; }
1392  if (ix >= inputValues->Nx) { ix = inputValues->Nx - 1; bflag = 0; }
1393  if (iy >= inputValues->Ny) { iy = inputValues->Ny - 1; bflag = 0; }
1394  if (iz < 0) { iz = 0; photon_status = DEAD; sleft = 0; }
1395  if (ix < 0) { ix = 0; bflag = 0; }
1396  if (iy < 0) { iy = 0; bflag = 0; }
1397  }
1398 
1399  // update pointer to tissue type
1400  i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy);
1401  } //(sv) /* same voxel */
1402  } while (sleft > 0); //do...while
1403 
1404  /**** SPIN
1405  Scatter photon into new trajectory defined by theta and psi.
1406  Theta is specified by cos(theta), which is determined
1407  based on the Henyey-Greenstein scattering function.
1408  Convert theta and psi into cosines ux, uy, uz.
1409  *****/
1410  /* Sample for costheta */
1411  while ((rnd = returnValue->RandomGen(1, 0, nullptr)) <= 0.0);
1412  if (inputValues->gVector[i] == 0.0)
1413  {
1414  costheta = 2.0 * rnd - 1.0;
1415  }
1416  else
1417  {
1418  double temp = (1.0 - inputValues->gVector[i] * inputValues->gVector[i])
1419  / (1.0 - inputValues->gVector[i] + 2 * inputValues->gVector[i] * rnd);
1420  costheta = (1.0 + inputValues->gVector[i] * inputValues->gVector[i] - temp*temp) / (2.0*inputValues->gVector[i]);
1421  }
1422  sintheta = sqrt(1.0 - costheta*costheta); /* sqrt() is faster than sin(). */
1423 
1424  /* Sample psi. */
1425  psi = 2.0*PI*returnValue->RandomGen(1, 0, nullptr);
1426  cospsi = cos(psi);
1427  if (psi < PI)
1428  sinpsi = sqrt(1.0 - cospsi*cospsi); /* sqrt() is faster than sin(). */
1429  else
1430  sinpsi = -sqrt(1.0 - cospsi*cospsi);
1431 
1432  /* New trajectory. */
1433  if (1 - fabs(uz) <= ONE_MINUS_COSZERO) { /* close to perpendicular. */
1434  uxx = sintheta * cospsi;
1435  uyy = sintheta * sinpsi;
1436  uzz = costheta * SIGN(uz); /* SIGN() is faster than division. */
1437  }
1438  else { /* usually use this option */
1439  temp = sqrt(1.0 - uz * uz);
1440  uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta;
1441  uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta;
1442  uzz = -sintheta * cospsi * temp + uz * costheta;
1443  }
1444 
1445  /* Update trajectory */
1446  ux = uxx;
1447  uy = uyy;
1448  uz = uzz;
1449 
1450  /**** CHECK ROULETTE
1451  If photon weight below THRESHOLD, then terminate photon using Roulette technique.
1452  Photon has CHANCE probability of having its weight increased by factor of 1/CHANCE,
1453  and 1-CHANCE probability of terminating.
1454  *****/
1455  if (W < THRESHOLD) {
1456  if (returnValue->RandomGen(1, 0, nullptr) <= CHANCE)
1457  W /= CHANCE;
1458  else photon_status = DEAD;
1459  }
1460  } while (photon_status == ALIVE); /* end STEP_CHECK_HOP_SPIN */
1461  /* if ALIVE, continue propagating */
1462  /* If photon DEAD, then launch new photon. */
1463  } while (photonIterator < photonsToSimulate); /* end RUN */
1464 
1465  returnValue->Nphotons += photonsToSimulate;
1466  } while (photonsToSimulate > 0);
1467 
1468  if (verbose) std::cout << "------------------------------------------------------" << std::endl;
1469  if (verbose) std::cout << "Thread " << thread << " is finished." << std::endl;
1470 }
#define ls
Definition: MitkMCxyz.cpp:57
int concurentThreadsSupported
Definition: MitkMCxyz.cpp:621
#define DEAD
Definition: MitkMCxyz.cpp:60
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define MZ
Definition: MitkMCxyz.cpp:287
std::string inputFilename
Definition: MitkMCxyz.cpp:625
#define ALIVE
Definition: MitkMCxyz.cpp:59
void setContributor(std::string contributor)
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
float yOffset
Definition: MitkMCxyz.cpp:622
void addArgument(const std::string &longarg, const std::string &shortarg, Type type, const std::string &argLabel, const std::string &argHelp=std::string(), const us::Any &defaultValue=us::Any(), bool optional=true, bool ignoreRest=false, bool deprecated=false, mitkCommandLineParser::Channel channel=mitkCommandLineParser::Channel::None)
float requestedSimulationTime
Definition: MitkMCxyz.cpp:620
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
std::string outputFilename
Definition: MitkMCxyz.cpp:626
static void info(const char *fmt,...)
Definition: svm.cpp:86
int main(int argc, char *argv[])
Definition: MitkMCxyz.cpp:630
mitk::pa::Probe::Pointer m_PhotoacousticProbe
Definition: MitkMCxyz.cpp:628
static IPropertyPersistence * GetPropertyPersistence(us::ModuleContext *context=us::GetModuleContext())
Get an IPropertyPersistence instance.
bool saveLegacy
Definition: MitkMCxyz.cpp:623
void runMonteCarlo(InputValues *inputValues, ReturnValues *returnValue, int thread, mitk::pa::MonteCarloThreadHandler::Pointer threadHandler)
Definition: MitkMCxyz.cpp:988
#define ONE_MINUS_COSZERO
Definition: MitkMCxyz.cpp:65
#define mitkThrow()
Definition: usAny.h:163
void setCategory(std::string category)
static Pointer New()
void setArgumentPrefix(const std::string &longPrefix, const std::string &shortPrefix)
std::string normalizationFilename
Definition: MitkMCxyz.cpp:624
bool verbose(false)
#define THRESHOLD
Definition: MitkMCxyz.cpp:61
int requestedNumberOfPhotons
Definition: MitkMCxyz.cpp:619
static void Save(const mitk::BaseData *data, const std::string &path, bool setPathProperty=false)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:774
bool interpretAsTime
Definition: MitkMCxyz.cpp:617
#define SIGN(x)
Definition: MitkMCxyz.cpp:64
#define SQR(x)
Definition: MitkMCxyz.cpp:63
MITKCORE_EXPORT const ScalarType eps
#define PI
Definition: MitkMCxyz.cpp:58
virtual bool AddInfo(const PropertyPersistenceInfo *info, bool overwrite=false)=0
Add persistence info for a specific base data property. If there is already a property info instance ...
void setTitle(std::string title)
int detector_x
Definition: MitkMCxyz.cpp:615
ImageReadAccessor class to get locked read access for a particular image part.
void setDescription(std::string description)
#define CHANCE
Definition: MitkMCxyz.cpp:62
void beginGroup(const std::string &description)
struct Location initLocation(int x, int y, int z, double absorb)
Definition: MitkMCxyz.cpp:75
#define FAC
Definition: MitkMCxyz.cpp:288
const void * GetData() const
Gives const access to the data.
#define MBIG
Definition: MitkMCxyz.cpp:285
int detector_z
Definition: MitkMCxyz.cpp:616
Class for defining the data type of pixels.
Definition: mitkPixelType.h:51
#define MSEED
Definition: MitkMCxyz.cpp:286
bool simulatePVFC
Definition: MitkMCxyz.cpp:618