Medical Imaging Interaction Toolkit  2018.4.99-36d69b77
Medical Imaging Interaction Toolkit
mitkPATissueGenerator.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 "mitkPATissueGenerator.h"
14 #include "mitkPAVector.h"
16 
17 mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData(
18  TissueGeneratorParameters::Pointer parameters)
19 {
20  MITK_DEBUG << "Initializing Empty Volume";
21 
22  const double RESAMPLING_FACTOR = 2;
23 
24  if (parameters->GetDoPartialVolume())
25  {
26  parameters->SetXDim(parameters->GetXDim() * RESAMPLING_FACTOR);
27  parameters->SetYDim(parameters->GetYDim() * RESAMPLING_FACTOR);
28  parameters->SetZDim(parameters->GetZDim() * RESAMPLING_FACTOR);
29  parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() * RESAMPLING_FACTOR);
30  parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() / RESAMPLING_FACTOR);
31  }
32 
33  std::mt19937 randomNumberGenerator;
34  std::random_device randomDevice;
35  if (parameters->GetUseRngSeed())
36  {
37  randomNumberGenerator.seed(parameters->GetRngSeed());
38  }
39  else
40  {
41  if (randomDevice.entropy() > 0.1)
42  {
43  randomNumberGenerator.seed(randomDevice());
44  }
45  else
46  {
47  randomNumberGenerator.seed(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
48  }
49  }
50 
51  auto generatedVolume = mitk::pa::InSilicoTissueVolume::New(parameters, &randomNumberGenerator);
52 
53  double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2;
54 
55  if(parameters->GetForceVesselsMoveAlongYDirection())
56  DIRECTION_VECTOR_INITIAL_VARIANCE = 0;
57 
58  std::uniform_int_distribution<int> randomNumVesselDistribution(
59  parameters->GetMinNumberOfVessels(), parameters->GetMaxNumberOfVessels());
60  std::uniform_real_distribution<double> randomBendingDistribution(
61  parameters->GetMinVesselBending(), parameters->GetMaxVesselBending());
62  std::uniform_real_distribution<double> randomAbsorptionDistribution(
63  parameters->GetMinVesselAbsorption(), parameters->GetMaxVesselAbsorption());
64  std::uniform_real_distribution<double> randomRadiusDistribution(
65  parameters->GetMinVesselRadiusInMillimeters(), parameters->GetMaxVesselRadiusInMillimeters());
66  std::uniform_real_distribution<double> randomScatteringDistribution(
67  parameters->GetMinVesselScattering(), parameters->GetMaxVesselScattering());
68  std::uniform_real_distribution<double> randomAnisotropyDistribution(
69  parameters->GetMinVesselAnisotropy(), parameters->GetMaxVesselAnisotropy());
70  std::uniform_int_distribution<int> borderTypeDistribution(0, 3);
71 
72  int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator);
73 
74  generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels);
75  generatedVolume->AddIntProperty("bifurcationFrequency", parameters->GetVesselBifurcationFrequency());
76 
77  MITK_INFO << "Simulating " << numberOfBloodVessels << " vessel structures";
78  for (int vesselNumber = 0; vesselNumber < numberOfBloodVessels; vesselNumber++)
79  {
80  Vector::Pointer initialPosition = Vector::New();
81  Vector::Pointer initialDirection = Vector::New();
82 
83  double initialRadius = randomRadiusDistribution(randomNumberGenerator) / parameters->GetVoxelSpacingInCentimeters() / 10;
84  std::stringstream radiusString;
85  radiusString << "vessel_" << vesselNumber + 1 << "_radius";
86  generatedVolume->AddDoubleProperty(radiusString.str(), initialRadius);
87 
88  double absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator);
89  std::stringstream absorptionString;
90  absorptionString << "vessel_" << vesselNumber + 1 << "_absorption";
91  generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient);
92 
93  double bendingFactor = randomBendingDistribution(randomNumberGenerator);
94  std::stringstream bendingString;
95  bendingString << "vessel_" << vesselNumber + 1 << "_bendingFactor";
96  generatedVolume->AddDoubleProperty(bendingString.str(), bendingFactor);
97 
98  double vesselScattering = randomScatteringDistribution(randomNumberGenerator);
99  std::stringstream scatteringString;
100  scatteringString << "vessel_" << vesselNumber + 1 << "_scattering";
101  generatedVolume->AddDoubleProperty(scatteringString.str(), vesselScattering);
102 
103  double vesselAnisotropy = randomAnisotropyDistribution(randomNumberGenerator);
104  std::stringstream anisotropyString;
105  anisotropyString << "vessel_" << vesselNumber + 1 << "_anisotropy";
106  generatedVolume->AddDoubleProperty(anisotropyString.str(), vesselAnisotropy);
107 
108  /*The vessel tree shall start at one of the 4 sides of the volume.
109  * The vessels will always be completely contained in the volume
110  * when starting to meander.
111  * They will meander in a direction perpendicular to the
112  * plane they started from (within the limits of the
113  * DIRECTION_VECTOR_INITIAL_VARIANCE)
114  */
115  int borderType = borderTypeDistribution(randomNumberGenerator);
116 
117  if(parameters->GetForceVesselsMoveAlongYDirection())
118  borderType = 2;
119 
120  switch (borderType)
121  {
122  case 0:
123  initialPosition->Randomize(0, 0, initialRadius, parameters->GetYDim() - initialRadius,
124  parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(),
125  parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator);
126  initialDirection->Randomize(1, 2, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE,
127  -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator);
128  break;
129  case 1:
130  initialPosition->Randomize(parameters->GetXDim() - 1, parameters->GetXDim() - 1, initialRadius, parameters->GetYDim() - initialRadius,
131  parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(),
132  parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator);
133  initialDirection->Randomize(-2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE,
134  -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator);
135  break;
136  case 2:
137  initialPosition->Randomize(initialRadius, parameters->GetXDim() - initialRadius, 0, 0,
138  parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(),
139  parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator);
140  initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, 1, 2,
141  -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator);
142  break;
143  case 3:
144  initialPosition->Randomize(initialRadius, parameters->GetXDim() - initialRadius, parameters->GetYDim() - 1, parameters->GetYDim() - 1,
145  parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(),
146  parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator);
147  initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -2, -1,
148  -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator);
149  break;
150  }
151 
152  initialDirection->Normalize();
153  MITK_INFO << initialPosition->GetElement(0) << " | " << initialPosition->GetElement(1) << " | " << initialPosition->GetElement(2);
154  MITK_INFO << initialDirection->GetElement(0) << " | " << initialDirection->GetElement(1) << " | " << initialDirection->GetElement(2);
155 
156  VesselProperties::Pointer vesselParams = VesselProperties::New();
157  vesselParams->SetDirectionVector(initialDirection);
158  vesselParams->SetPositionVector(initialPosition);
159  vesselParams->SetRadiusInVoxel(initialRadius);
160  vesselParams->SetAbsorptionCoefficient(absorptionCoefficient);
161  vesselParams->SetScatteringCoefficient(vesselScattering);
162  vesselParams->SetAnisotopyCoefficient(vesselAnisotropy);
163  vesselParams->SetBifurcationFrequency(parameters->GetVesselBifurcationFrequency());
164  vesselParams->SetDoPartialVolume(parameters->GetDoPartialVolume());
165 
166  VesselTree::Pointer vesselTree = VesselTree::New(vesselParams);
167 
168  while (!vesselTree->IsFinished())
169  {
170  vesselTree->Step(generatedVolume, parameters->GetCalculateNewVesselPositionCallback(), bendingFactor, &randomNumberGenerator);
171  }
172  }
173 
174  if (parameters->GetDoPartialVolume())
175  {
176  VolumeManipulator::RescaleImage(generatedVolume, (1.0 / RESAMPLING_FACTOR));
177  parameters->SetXDim(parameters->GetXDim() / RESAMPLING_FACTOR);
178  parameters->SetYDim(parameters->GetYDim() / RESAMPLING_FACTOR);
179  parameters->SetZDim(parameters->GetZDim() / RESAMPLING_FACTOR);
180  parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() / RESAMPLING_FACTOR);
181  parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() * RESAMPLING_FACTOR);
182  }
183 
184  generatedVolume->FinalizeVolume();
185 
186  return generatedVolume;
187 }
188 
189 mitk::pa::InSilicoTissueGenerator::InSilicoTissueGenerator()
190 {
191 }
192 
193 mitk::pa::InSilicoTissueGenerator::~InSilicoTissueGenerator()
194 {
195 }
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define MITK_DEBUG
Definition: mitkLogMacros.h:22
static InSilicoTissueVolume::Pointer New(mitk::pa::Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList)
static InSilicoTissueVolume::Pointer GenerateInSilicoData(TissueGeneratorParameters::Pointer parameters)
GenerateInSilicoData This method will return a InSilicoTissueVolume created in terms of the given par...
static void RescaleImage(InSilicoTissueVolume::Pointer image, double ratio)
static Pointer New()