Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkPAInSilicoTissueVolume.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 
14 #include <mitkProperties.h>
15 #include <mitkCoreServices.h>
17 #include <itkImage.h>
18 #include <itkDiscreteGaussianImageFilter.h>
19 #include <mitkImageCast.h>
20 #include <mitkImageToItk.h>
21 #include <chrono>
22 
23 mitk::pa::InSilicoTissueVolume::InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters, std::mt19937* rng)
24 {
25  {
26  unsigned int xDim = parameters->GetXDim();
27  unsigned int yDim = parameters->GetYDim();
28  unsigned int zDim = parameters->GetZDim();
29  m_TDim = 4;
30  unsigned int size = xDim * yDim * zDim;
31  auto* absorptionArray = new double[size];
32  auto* scatteringArray = new double[size];
33  auto* anisotropyArray = new double[size];
34  auto* segmentationArray = new double[size];
35 
36  m_InitialBackgroundAbsorption = (parameters->GetMinBackgroundAbsorption() + parameters->GetMaxBackgroundAbsorption()) / 2;
37  m_Rng = rng;
38 
39  for (unsigned int index = 0; index < size; index++)
40  {
41  absorptionArray[index] = m_InitialBackgroundAbsorption;
42  scatteringArray[index] = parameters->GetBackgroundScattering();
43  anisotropyArray[index] = parameters->GetBackgroundAnisotropy();
44  segmentationArray[index] = SegmentationType::BACKGROUND;
45  }
46 
47  m_AbsorptionVolume = Volume::New(absorptionArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters());
48  m_ScatteringVolume = Volume::New(scatteringArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters());
49  m_AnisotropyVolume = Volume::New(anisotropyArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters());
50  m_SegmentationVolume = Volume::New(segmentationArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters());
51  }
52 
53  m_TissueParameters = parameters;
55  UpdatePropertyList();
56 }
57 
58 void mitk::pa::InSilicoTissueVolume::UpdatePropertyList()
59 {
60  //Set properties
61  AddIntProperty("mcflag", m_TissueParameters->GetMCflag());
62  AddIntProperty("launchflag", m_TissueParameters->GetMCLaunchflag());
63  AddIntProperty("boundaryflag", m_TissueParameters->GetMCBoundaryflag());
64  AddDoubleProperty("launchPointX", m_TissueParameters->GetMCLaunchPointX());
65  AddDoubleProperty("launchPointY", m_TissueParameters->GetMCLaunchPointY());
66  AddDoubleProperty("launchPointZ", m_TissueParameters->GetMCLaunchPointZ());
67  AddDoubleProperty("focusPointX", m_TissueParameters->GetMCFocusPointX());
68  AddDoubleProperty("focusPointY", m_TissueParameters->GetMCFocusPointY());
69  AddDoubleProperty("focusPointZ", m_TissueParameters->GetMCFocusPointZ());
70  AddDoubleProperty("trajectoryVectorX", m_TissueParameters->GetMCTrajectoryVectorX());
71  AddDoubleProperty("trajectoryVectorY", m_TissueParameters->GetMCTrajectoryVectorY());
72  AddDoubleProperty("trajectoryVectorZ", m_TissueParameters->GetMCTrajectoryVectorZ());
73  AddDoubleProperty("radius", m_TissueParameters->GetMCRadius());
74  AddDoubleProperty("waist", m_TissueParameters->GetMCWaist());
75  AddDoubleProperty("partialVolume", m_TissueParameters->GetDoPartialVolume());
76  AddDoubleProperty("standardTissueAbsorptionMin", m_TissueParameters->GetMinBackgroundAbsorption());
77  AddDoubleProperty("standardTissueAbsorptionMax", m_TissueParameters->GetMaxBackgroundAbsorption());
78  AddDoubleProperty("standardTissueScattering", m_TissueParameters->GetBackgroundScattering());
79  AddDoubleProperty("standardTissueAnisotropy", m_TissueParameters->GetBackgroundAnisotropy());
80  AddDoubleProperty("airThickness", m_TissueParameters->GetAirThicknessInMillimeters());
81  AddDoubleProperty("skinThickness", m_TissueParameters->GetSkinThicknessInMillimeters());
82 }
83 
85  Volume::Pointer absorptionVolume,
86  Volume::Pointer scatteringVolume,
87  Volume::Pointer anisotropyVolume,
88  Volume::Pointer segmentationVolume,
89  TissueGeneratorParameters::Pointer tissueParameters,
90  mitk::PropertyList::Pointer propertyList)
91 {
92  m_AbsorptionVolume = absorptionVolume;
93  m_ScatteringVolume = scatteringVolume;
94  m_AnisotropyVolume = anisotropyVolume;
95  m_SegmentationVolume = segmentationVolume;
96  m_TissueParameters = tissueParameters;
97  m_PropertyList = propertyList;
98  if (m_SegmentationVolume.IsNotNull())
99  m_TDim = 4;
100  else
101  m_TDim = 3;
102 }
103 
105 {
106  return m_AbsorptionVolume->GetSpacing();
107 }
108 
110 {
111  m_AbsorptionVolume->SetSpacing(spacing);
112  m_ScatteringVolume->SetSpacing(spacing);
113  m_AnisotropyVolume->SetSpacing(spacing);
114  if (m_SegmentationVolume.IsNotNull())
115  {
116  m_SegmentationVolume->SetSpacing(spacing);
117  }
118 }
119 
120 void mitk::pa::InSilicoTissueVolume::AddDoubleProperty(std::string label, double value)
121 {
122  m_PropertyList->SetDoubleProperty(label.c_str(), value);
124 }
125 
126 void mitk::pa::InSilicoTissueVolume::AddIntProperty(std::string label, int value)
127 {
128  m_PropertyList->SetIntProperty(label.c_str(), value);
130 }
131 
133 {
134  mitk::Image::Pointer resultImage = mitk::Image::New();
135  mitk::PixelType TPixel = mitk::MakeScalarPixelType<double>();
136  auto* dimensionsOfImage = new unsigned int[4];
137 
138  // Copy dimensions
139  dimensionsOfImage[0] = m_TissueParameters->GetYDim();
140  dimensionsOfImage[1] = m_TissueParameters->GetXDim();
141  dimensionsOfImage[2] = m_TissueParameters->GetZDim();
142  dimensionsOfImage[3] = m_TDim;
143 
144  resultImage->Initialize(TPixel, 4, dimensionsOfImage, 1);
145 
146  mitk::Vector3D spacing;
147  spacing.Fill(m_TissueParameters->GetVoxelSpacingInCentimeters());
148  resultImage->SetSpacing(spacing);
149 
150  resultImage->SetImportVolume(m_AbsorptionVolume->GetData(), 0, 0, mitk::Image::CopyMemory);
151  resultImage->SetImportVolume(m_ScatteringVolume->GetData(), 1, 0, mitk::Image::CopyMemory);
152  resultImage->SetImportVolume(m_AnisotropyVolume->GetData(), 2, 0, mitk::Image::CopyMemory);
153  if (m_TDim > 3)
154  {
155  resultImage->SetImportVolume(m_SegmentationVolume->GetData(), 3, 0, mitk::Image::CopyMemory);
156  }
157 
158  resultImage->SetPropertyList(m_PropertyList);
159 
160  return resultImage;
161 }
162 
163 mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueVolume::New(
164  Volume::Pointer absorptionVolume,
165  Volume::Pointer scatteringVolume,
166  Volume::Pointer anisotropyVolume,
167  Volume::Pointer segmentationVolume,
168  TissueGeneratorParameters::Pointer tissueParameters,
169  mitk::PropertyList::Pointer propertyList)
170 {
171  InSilicoTissueVolume::Pointer smartPtr = new InSilicoTissueVolume(
172  absorptionVolume, scatteringVolume, anisotropyVolume, segmentationVolume,
173  tissueParameters, propertyList);
174  smartPtr->UnRegister();
175  return smartPtr;
176 }
177 
179 {
180  m_AbsorptionVolume = nullptr;
181  m_ScatteringVolume = nullptr;
182  m_AnisotropyVolume = nullptr;
183  m_SegmentationVolume = nullptr;
184  m_TissueParameters = nullptr;
185  m_PropertyList = nullptr;
186 }
187 
188 void mitk::pa::InSilicoTissueVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy)
189 {
190  if (IsInsideVolume(x, y, z))
191  {
192  m_AbsorptionVolume->SetData(absorption, x, y, z);
193  m_ScatteringVolume->SetData(scattering, x, y, z);
194  m_AnisotropyVolume->SetData(anisotropy, x, y, z);
195  }
196 }
197 
198 void mitk::pa::InSilicoTissueVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType)
199 {
200  if (IsInsideVolume(x, y, z))
201  {
202  m_AbsorptionVolume->SetData(absorption, x, y, z);
203  m_ScatteringVolume->SetData(scattering, x, y, z);
204  m_AnisotropyVolume->SetData(anisotropy, x, y, z);
205  if (m_SegmentationVolume.IsNotNull())
206  {
207  m_SegmentationVolume->SetData(segmentType, x, y, z);
208  }
209  }
210 }
211 
213 {
214  return x >= 0 && x < m_TissueParameters->GetXDim() &&
215  y >= 0 && y < m_TissueParameters->GetYDim() &&
216  z >= 0 && z < m_TissueParameters->GetZDim();
217 }
218 
220 {
221  return m_AbsorptionVolume;
222 }
223 
225 {
226  return m_SegmentationVolume;
227 }
228 
230 {
231  AddSkinAndAirLayers();
232 
233  // If specified, randomize all tissue parameters
234  if (m_TissueParameters->GetRandomizePhysicalProperties())
235  {
237  m_TissueParameters->GetRngSeed(),
238  m_TissueParameters->GetRandomizePhysicalPropertiesPercentage());
239  }
240 
241  unsigned int xDim = m_TissueParameters->GetXDim();
242  unsigned int yDim = m_TissueParameters->GetYDim();
243  unsigned int zDim = m_TissueParameters->GetZDim();
244 
245  std::uniform_real_distribution<double> randomBackgroundAbsorptionDistribution(
246  m_TissueParameters->GetMinBackgroundAbsorption(), m_TissueParameters->GetMaxBackgroundAbsorption());
247 
248  for (unsigned int z = 0; z < zDim; z++)
249  {
250  for (unsigned int y = 0; y < yDim; y++)
251  {
252  for (unsigned int x = 0; x < xDim; x++)
253  {
254  if (fabs(m_AbsorptionVolume->GetData(x, y, z) - m_InitialBackgroundAbsorption) < mitk::eps)
255  {
256  m_AbsorptionVolume->SetData(randomBackgroundAbsorptionDistribution(*m_Rng), x, y, z);
257  }
258  }
259  }
260  }
261 }
262 
263 void mitk::pa::InSilicoTissueVolume::AddSkinAndAirLayers()
264 {
265  //Calculate the index location according to thickness in cm
266  double airvoxel = (m_TissueParameters->GetAirThicknessInMillimeters() / m_TissueParameters->GetVoxelSpacingInCentimeters()) / 10;
267  double skinvoxel = airvoxel + (m_TissueParameters->GetSkinThicknessInMillimeters() / m_TissueParameters->GetVoxelSpacingInCentimeters()) / 10;
268 
269  for (int y = 0; y < m_TissueParameters->GetYDim(); y++)
270  {
271  for (int x = 0; x < m_TissueParameters->GetXDim(); x++)
272  {
273  // Add air from index 0 to airvoxel
274  if (m_TissueParameters->GetAirThicknessInMillimeters() > mitk::eps)
275  {
276  FillZLayer(x, y, 0, airvoxel,
277  m_TissueParameters->GetAirAbsorption(),
278  m_TissueParameters->GetAirScattering(),
279  m_TissueParameters->GetAirAnisotropy(),
280  SegmentationType::AIR);
281  }
282 
283  //Add skin from index airvoxel to skinvoxel
284  if (m_TissueParameters->GetSkinThicknessInMillimeters() > mitk::eps)
285  {
286  FillZLayer(x, y, airvoxel, skinvoxel,
287  m_TissueParameters->GetSkinAbsorption(),
288  m_TissueParameters->GetSkinScattering(),
289  m_TissueParameters->GetSkinAnisotropy(),
290  SegmentationType::SKIN);
291  }
292  }
293  }
294 }
295 
296 void mitk::pa::InSilicoTissueVolume::FillZLayer(int x, int y, double startIdx, double endIdx,
297  double absorption, double scattering, double anisotropy,
298  SegmentationType segmentationType)
299 {
300  for (int z = startIdx; z < endIdx; z++)
301  {
302  if (IsInsideVolume(x, y, z))
303  {
304  if (endIdx - z < 1)
305  {
306  //Simulate partial volume effects
307  m_AbsorptionVolume->SetData((1 - (endIdx - z)) *
308  m_AbsorptionVolume->GetData(x, y, z) + (endIdx - z) * absorption, x, y, z);
309  m_ScatteringVolume->SetData((1 - (endIdx - z)) *
310  m_ScatteringVolume->GetData(x, y, z) + (endIdx - z) * scattering, x, y, z);
311  m_AnisotropyVolume->SetData((1 - (endIdx - z)) *
312  m_AnisotropyVolume->GetData(x, y, z) + (endIdx - z) * anisotropy, x, y, z);
313  if (endIdx - z > 0.5)
314  {
315  //Only put the segmentation label if more than half of the partial volume is the wanted tissue type
316  if (m_SegmentationVolume.IsNotNull())
317  {
318  m_SegmentationVolume->SetData(segmentationType, x, y, z);
319  }
320  }
321  }
322  else
323  {
324  m_AbsorptionVolume->SetData(absorption, x, y, z);
325  m_ScatteringVolume->SetData(scattering, x, y, z);
326  m_AnisotropyVolume->SetData(anisotropy, x, y, z);
327  if (m_SegmentationVolume.IsNotNull())
328  {
329  m_SegmentationVolume->SetData(segmentationType, x, y, z);
330  }
331  }
332  }
333  }
334 }
335 
336 void mitk::pa::InSilicoTissueVolume::RandomizeTissueCoefficients(long rngSeed, bool useRngSeed, double percentage)
337 {
338  std::mt19937 rng;
339  std::random_device randomDevice;
340  if (useRngSeed)
341  {
342  rng.seed(rngSeed);
343  }
344  else
345  {
346  if (randomDevice.entropy() > 0.1)
347  {
348  rng.seed(randomDevice());
349  }
350  else
351  {
352  rng.seed(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
353  }
354  }
355  std::normal_distribution<> percentageDistribution(1, percentage / 100);
356 
357  for (int y = 0; y < m_TissueParameters->GetYDim(); y++)
358  {
359  for (int x = 0; x < m_TissueParameters->GetXDim(); x++)
360  {
361  for (int z = 0; z < m_TissueParameters->GetZDim(); z++)
362  {
363  m_AbsorptionVolume->SetData(m_AbsorptionVolume->GetData(x, y, z) * percentageDistribution(rng), x, y, z);
364  m_ScatteringVolume->SetData(m_ScatteringVolume->GetData(x, y, z) * percentageDistribution(rng), x, y, z);
365  }
366  }
367  }
368 }
369 
371 {
372  return m_ScatteringVolume;
373 }
374 
376 {
377  return m_AnisotropyVolume;
378 }
379 
381 {
382  m_AbsorptionVolume = volume;
383 }
384 
386 {
387  m_ScatteringVolume = volume;
388 }
389 
391 {
392  m_AnisotropyVolume = volume;
393 }
394 
396 {
397  m_SegmentationVolume = volume;
398 }
void SetAbsorptionVolume(Volume::Pointer volume)
mitk::pa::Volume::Pointer m_ScatteringVolume
static Pointer New()
void RandomizeTissueCoefficients(long rngSeed, bool useRngSeed, double percentage)
mitk::pa::Volume::Pointer m_AnisotropyVolume
TissueGeneratorParameters::Pointer m_TissueParameters
bool IsInsideVolume(int x, int y, int z)
IsInsideVolume.
void SetSegmentationVolume(Volume::Pointer volume)
mitk::PropertyList::Pointer m_PropertyList
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)
mitk::Image::Pointer ConvertToMitkImage()
ConvertToMitkImage.
void SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType)
SetVolumeValues sets the values for aborption, scattering and anisotropy at the specified voxel locat...
static IPropertyPersistence * GetPropertyPersistence(us::ModuleContext *context=us::GetModuleContext())
Get an IPropertyPersistence instance.
void SetAnisotropyVolume(Volume::Pointer volume)
static Pointer New()
void SetScatteringVolume(Volume::Pointer volume)
void AddDoubleProperty(std::string label, double value)
AddDoubleProperty adds a persistent property to the volume, which will be exported to the mitk image...
MITKCORE_EXPORT const ScalarType eps
mitk::pa::Volume::Pointer m_SegmentationVolume
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 ...
InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters, std::mt19937 *rng)
static Volume::Pointer New(double *data, unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing)
returns smartpointer reference to a new instance of this objects. The given data array will be freed ...
void AddIntProperty(std::string label, int value)
AddIntProperty adds a persistent property to the volume, which will be exported to the mitk image...
Class for defining the data type of pixels.
Definition: mitkPixelType.h:51
mitk::pa::Volume::Pointer m_AbsorptionVolume