Medical Imaging Interaction Toolkit  2018.4.99-9a29ffc6
Medical Imaging Interaction Toolkit
mitkPALightSource.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 
13 #include "mitkPALightSource.h"
14 #include <cmath>
15 
17  m_IsValid(false)
18 {
19 }
20 
21 mitk::pa::LightSource::LightSource(TiXmlElement* element, bool verbose) :
22  m_IsValid(true),
23  m_Verbose(verbose)
24 {
25  ParseEnergy(element);
26  ParsePhotonSpawnArea(element);
27  ParsePhotonDirection(element);
28 
29  if (m_IsValid)
30  {
31  if (m_Verbose)
32  std::cout << "Successfully created LightSource" << std::endl;
33  }
34  else
35  {
36  if (m_Verbose)
37  std::cout << "Failed creating LightSource." << std::endl;
38  }
39 }
40 
42 {
43 }
44 
46 {
47  TransformResult result;
48  result.z0 = sqrt(-2.0 * log(u1)) * cos(TWO_PI * u2) * sigma + mu;
49  result.z1 = sqrt(-2.0 * log(u1)) * sin(TWO_PI * u2) * sigma + mu;
50  return result;
51 }
52 
54 {
55  TiXmlElement* direction = element->FirstChildElement(XML_TAG_PHOTON_DIRECTION);
56  if (direction)
57  {
58  ParseAngle(direction, XML_TAG_X_ANGLE);
59  }
60  else
61  {
62  if (m_Verbose)
63  std::cerr << "No \"" << XML_TAG_X_ANGLE << "\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl;
64 
65  m_AngleXMinimum = 0;
66  m_AngleXMaximum = 0;
67  m_AngleXMode = DistributionMode::UNIFORM;
68  }
69 
70  direction = element->FirstChildElement(XML_TAG_PHOTON_DIRECTION);
71  if (direction)
72  {
73  ParseAngle(direction, XML_TAG_Y_ANGLE);
74  }
75  else
76  {
77  if (m_Verbose)
78  std::cerr << "No \"" << XML_TAG_Y_ANGLE << "\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl;
79 
80  m_AngleYMinimum = 0;
81  m_AngleYMaximum = 0;
82  m_AngleYMode = DistributionMode::UNIFORM;
83  }
84 }
85 
86 void mitk::pa::LightSource::ParseAngle(TiXmlElement* direction, std::string angle)
87 {
88  double minimum;
89  double maximum;
90  DistributionMode mode = DistributionMode::GAUSSIAN;
91 
92  if (m_Verbose)
93  std::cout << "Parsing " << angle << std::endl;
94  TiXmlElement* angleElement = direction->FirstChildElement(angle);
95  if (angleElement)
96  {
97  TiXmlElement* angleMin = angleElement->FirstChildElement(XML_TAG_MINIMUM);
98  if (angleMin)
99  {
100  std::string angleMinText = angleMin->GetText();
101  minimum = std::stod(angleMinText);
102  if (m_Verbose)
103  std::cout << "Setting min=" << minimum << std::endl;
104  }
105  else
106  {
107  if (m_Verbose)
108  std::cerr << "No \"" << XML_TAG_MINIMUM << "\" tag in xml. Setting min=0" << std::endl;
109  minimum = 0;
110  }
111 
112  TiXmlElement* angleMax = angleElement->FirstChildElement(XML_TAG_MAXIMUM);
113  if (angleMax)
114  {
115  std::string angleMaxText = angleMax->GetText();
116  maximum = std::stod(angleMaxText);
117  if (m_Verbose)
118  std::cout << "Setting max=" << maximum << std::endl;
119  }
120  else
121  {
122  if (m_Verbose)
123  std::cerr << "No \"" << XML_TAG_MAXIMUM << "\" tag in xml. Setting max=0" << std::endl;
124  maximum = 0;
125  }
126 
127  TiXmlElement* angleMode = angleElement->FirstChildElement(XML_TAG_MODE);
128  if (angleMode)
129  {
130  std::string angleModeText = angleMode->GetText();
131  if (strcmp("UNIFORM", angleModeText.c_str()) == 0)
132  {
133  mode = DistributionMode::UNIFORM;
134  if (m_Verbose)
135  std::cout << "Setting mode=UNIFORM" << std::endl;
136  }
137  else if (strcmp("GAUSSIAN", angleModeText.c_str()) == 0)
138  {
139  mode = DistributionMode::GAUSSIAN;
140  if (m_Verbose)
141  std::cout << "Setting mode=GAUSSIAN" << std::endl;
142  }
143  }
144  else
145  {
146  if (m_Verbose)
147  std::cerr << "No \"" << XML_TAG_MODE << "\" tag in xml. Setting mode=UNIFORM" << std::endl;
148  mode = DistributionMode::UNIFORM;
149  }
150  }
151  else
152  {
153  if (m_Verbose)
154  std::cerr << "No \"" << angle << "\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl;
155 
156  maximum = 0;
157  minimum = 0;
158  mode = DistributionMode::UNIFORM;
159  }
160 
161  if (strcmp(XML_TAG_X_ANGLE.c_str(), angle.c_str()) == 0)
162  {
163  m_AngleXMinimum = minimum;
164  m_AngleXMaximum = maximum;
165  m_AngleXMode = mode;
166  }
167  else if (strcmp(XML_TAG_Y_ANGLE.c_str(), angle.c_str()) == 0)
168  {
169  m_AngleYMinimum = minimum;
170  m_AngleYMaximum = maximum;
171  m_AngleYMode = mode;
172  }
173 }
174 
175 void mitk::pa::LightSource::ParseEnergy(TiXmlElement* element)
176 {
177  TiXmlElement* energy = element->FirstChildElement(XML_TAG_ENERGY);
178  if (energy)
179  {
180  std::string energyText = energy->GetText();
181  m_Energy = std::stod(energyText);
182  if (m_Verbose)
183  std::cout << "Setting energy=" << m_Energy;
184  }
185  else
186  {
187  if (m_Verbose)
188  std::cerr << "No \"" << XML_TAG_ENERGY << "\" field in xml. Setting Energy=1" << std::endl;
189  m_Energy = 1.0;
190  }
191 }
192 
194 {
195  TiXmlElement* spawnArea = element->FirstChildElement("PhotonSpawnArea");
196  if (spawnArea)
197  {
198  TiXmlElement* spawnType = spawnArea->FirstChildElement(XML_TAG_SPAWN_TYPE);
199  if (spawnType)
200  {
201  std::string spawnTypeText = spawnType->GetText();
202  if (strcmp(XML_TAG_SPAWN_TYPE_POINT.c_str(), spawnTypeText.c_str()) == 0)
203  {
204  m_SpawnType = SpawnType::POINT;
205  if (m_Verbose)
206  std::cout << "Setting " << XML_TAG_SPAWN_TYPE << " = " << XML_TAG_SPAWN_TYPE_POINT << std::endl;
207  }
208  else if (strcmp(XML_TAG_SPAWN_TYPE_RECTANGLE.c_str(), spawnTypeText.c_str()) == 0)
209  {
210  m_SpawnType = SpawnType::RECTANGLE;
211  if (m_Verbose)
212  std::cout << "Setting " << XML_TAG_SPAWN_TYPE << " = " << XML_TAG_SPAWN_TYPE_RECTANGLE << std::endl;
213  }
214  else if (strcmp(XML_TAG_SPAWN_TYPE_CIRCLE.c_str(), spawnTypeText.c_str()) == 0)
215  {
216  m_SpawnType = SpawnType::CIRCLE;
217  if (m_Verbose)
218  std::cout << "Setting " << XML_TAG_SPAWN_TYPE << " = " << XML_TAG_SPAWN_TYPE_CIRCLE << std::endl;
219  }
220  else
221  {
222  std::cerr << "The provided SpawnType (" << spawnTypeText << ") did not match any available spawn type. Light source is not valid." << std::endl;
223  m_IsValid = false;
224  }
225  }
226  else
227  {
228  std::cerr << "The \"" << XML_TAG_SPAWN_TYPE << "\" element was not provided for this light source. Light source is not valid." << std::endl;
229  m_IsValid = false;
230  }
231 
232  TiXmlElement* xLocation = spawnArea->FirstChildElement(XML_TAG_X);
233  if (xLocation)
234  {
235  std::string xLocationText = xLocation->GetText();
236  m_SpawnLocationX = std::stod(xLocationText);
237  if (m_Verbose)
238  std::cout << "Setting " << XML_TAG_X << "=" << m_SpawnLocationX;
239  }
240  else
241  {
242  if (m_Verbose)
243  std::cerr << "No \"" << XML_TAG_X << "\" field in xml. Setting " << XML_TAG_X << "=0" << std::endl;
244  m_SpawnLocationX = 0;
245  }
246 
247  TiXmlElement* yLocation = spawnArea->FirstChildElement(XML_TAG_Y);
248  if (yLocation)
249  {
250  std::string yLocationText = yLocation->GetText();
251  m_SpawnLocationY = std::stod(yLocationText);
252  if (m_Verbose)
253  std::cout << "Setting " << XML_TAG_Y << "=" << m_SpawnLocationY;
254  }
255  else
256  {
257  if (m_Verbose)
258  std::cerr << "No \"" << XML_TAG_Y << "\" field in xml. Setting " << XML_TAG_Y << "=0" << std::endl;
259  m_SpawnLocationY = 0;
260  }
261 
262  TiXmlElement* zLocation = spawnArea->FirstChildElement(XML_TAG_Z);
263  if (zLocation)
264  {
265  std::string zLocationText = zLocation->GetText();
266  m_SpawnLocationZ = std::stod(zLocationText);
267  if (m_Verbose)
268  std::cout << "Setting " << XML_TAG_Z << "=" << m_SpawnLocationZ;
269  }
270  else
271  {
272  if (m_Verbose)
273  std::cerr << "No \"" << XML_TAG_Z << "\" field in xml. Setting " << XML_TAG_Z << "=0.1" << std::endl;
274  m_SpawnLocationZ = 0.1;
275  }
276 
277  TiXmlElement* rLocation = spawnArea->FirstChildElement(XML_TAG_R);
278  if (rLocation)
279  {
280  std::string rLocationText = rLocation->GetText();
281  m_SpawnLocationRadius = std::stod(rLocationText);
282  if (m_Verbose)
283  std::cout << "Setting " << XML_TAG_R << "=" << m_SpawnLocationRadius;
284  }
285  else
286  {
287  if (m_Verbose)
288  std::cerr << "No \"" << XML_TAG_R << "\" field in xml. Setting " << XML_TAG_R << "=0" << std::endl;
290  }
291 
292  TiXmlElement* xLength = spawnArea->FirstChildElement(XML_TAG_X_LENGTH);
293  if (xLength)
294  {
295  std::string xLengthText = xLength->GetText();
296  m_SpawnLocationXLength = std::stod(xLengthText);
297  if (m_Verbose)
298  std::cout << "Setting " << XML_TAG_X_LENGTH << "=" << m_SpawnLocationXLength << std::endl;
299  }
300  else
301  {
302  if (m_Verbose)
303  std::cerr << "No \"" << XML_TAG_X_LENGTH << "\" field in xml. Setting " << XML_TAG_X_LENGTH << "=0" << std::endl;
305  }
306 
307  TiXmlElement* yLength = spawnArea->FirstChildElement(XML_TAG_Y_LENGTH);
308  if (yLength)
309  {
310  std::string yLengthText = yLength->GetText();
311  m_SpawnLocationYLength = std::stod(yLengthText);
312  if (m_Verbose)
313  std::cout << "Setting " << XML_TAG_Y_LENGTH << "=" << m_SpawnLocationYLength << std::endl;
314  }
315  else
316  {
317  if (m_Verbose)
318  std::cerr << "No \"" << XML_TAG_Y_LENGTH << "\" field in xml. Setting " << XML_TAG_Y_LENGTH << "=0" << std::endl;
320  }
321 
322  TiXmlElement* zLength = spawnArea->FirstChildElement(XML_TAG_Z_LENGTH);
323  if (zLength)
324  {
325  std::string zLengthText = zLength->GetText();
326  m_SpawnLocationZLength = std::stod(zLengthText);
327  if (m_Verbose)
328  std::cout << "Setting " << XML_TAG_Z_LENGTH << "=" << m_SpawnLocationZLength << std::endl;
329  }
330  else
331  {
332  if (m_Verbose)
333  std::cerr << "No \"" << XML_TAG_Z_LENGTH << "\" field in xml. Setting " << XML_TAG_Z_LENGTH << "=0" << std::endl;
335  }
336  }
337  else
338  m_IsValid = false;
339 }
340 
342  double rnd4, double rnd5, double gau1, double gau2)
343 {
344  PhotonInformation returnValue;
345 
346  switch (m_SpawnType)
347  {
348  case POINT:
349  returnValue.xPosition = m_SpawnLocationX;
350  returnValue.yPosition = m_SpawnLocationY;
351  returnValue.zPosition = m_SpawnLocationZ;
352  break;
353  case RECTANGLE:
354  returnValue.xPosition = m_SpawnLocationX + rnd3 * m_SpawnLocationXLength;
355  returnValue.yPosition = m_SpawnLocationY + rnd4 * m_SpawnLocationYLength;
356  returnValue.zPosition = m_SpawnLocationZ + rnd5 * m_SpawnLocationZLength;
357  break;
358  case CIRCLE:
359  double radius = rnd3 * m_SpawnLocationRadius;
360  double angle = rnd4 * TWO_PI;
361 
362  returnValue.xPosition = m_SpawnLocationX + radius * cos(angle);
363  returnValue.yPosition = m_SpawnLocationY + radius * sin(angle);
364  returnValue.zPosition = m_SpawnLocationZ;
365  break;
366  }
367 
368  switch (m_AngleXMode)
369  {
370  case UNIFORM:
371  returnValue.xAngle = rnd1 * (m_AngleXMaximum - m_AngleXMinimum) + m_AngleXMinimum;
372  break;
373  case GAUSSIAN:
375  returnValue.xAngle = trans.z0;
376  break;
377  }
378 
379  switch (m_AngleYMode)
380  {
381  case UNIFORM:
382  returnValue.yAngle = rnd2 * (m_AngleYMaximum - m_AngleYMinimum) + m_AngleYMinimum;
383  break;
384  case GAUSSIAN:
386  returnValue.yAngle = trans.z1;
387  break;
388  }
389 
390  if ((returnValue.xAngle*returnValue.xAngle + returnValue.yAngle*returnValue.yAngle) > 1)
391  {
392  double unify = sqrt(returnValue.xAngle*returnValue.xAngle + returnValue.yAngle*returnValue.yAngle)*1.001;
393  returnValue.xAngle = returnValue.xAngle / unify;
394  returnValue.yAngle = returnValue.yAngle / unify;
395  }
396 
397  returnValue.zAngle = sqrt(1 - returnValue.xAngle*returnValue.xAngle - returnValue.yAngle*returnValue.yAngle);
398 
399  if (m_Verbose)
400  std::cout << "Created a new photon at (" << returnValue.xPosition << "|" << returnValue.yPosition << "|" <<
401  returnValue.zPosition << ") with angle (" << returnValue.xAngle << "|" << returnValue.yAngle << "|" <<
402  returnValue.zAngle << ")" << std::endl;
403 
404  return returnValue;
405 }
406 
408 {
409  return m_IsValid;
410 }
void ParseEnergy(TiXmlElement *element)
const std::string XML_TAG_SPAWN_TYPE
const std::string XML_TAG_MAXIMUM
const std::string XML_TAG_ENERGY
const std::string XML_TAG_R
TransformResult BoxMuellerTransform(double u1, double u2, double mu, double sigma)
void ParsePhotonDirection(TiXmlElement *element)
const std::string XML_TAG_SPAWN_TYPE_POINT
DistributionMode m_AngleXMode
const std::string XML_TAG_Y_LENGTH
void ParsePhotonSpawnArea(TiXmlElement *element)
bool verbose(false)
const std::string XML_TAG_SPAWN_TYPE_CIRCLE
const std::string XML_TAG_X_LENGTH
const std::string XML_TAG_Z_LENGTH
void ParseAngle(TiXmlElement *direction, std::string angle)
const std::string XML_TAG_X
const std::string XML_TAG_PHOTON_DIRECTION
const std::string XML_TAG_MODE
const std::string XML_TAG_X_ANGLE
DistributionMode m_AngleYMode
const std::string XML_TAG_Y_ANGLE
const std::string XML_TAG_Z
const std::string XML_TAG_Y
PhotonInformation GetNextPhoton(double rnd1, double rnd2, double rnd3, double rnd4, double rnd5, double gau1, double gau2)
const std::string XML_TAG_SPAWN_TYPE_RECTANGLE
const std::string XML_TAG_MINIMUM