Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
PABeamformingTool.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 <mitkCommon.h>
14 #include <chrono>
15 #include <mitkIOUtil.h>
16 #include <mitkCommandLineParser.h>
17 #include <mitkException.h>
18 
22 
23 #include <itksys/SystemTools.hxx>
24 #include <tinyxml/tinyxml.h>
25 
26 struct InputParameters
27 {
28  mitk::Image::Pointer inputImage;
29  std::string outputFilename;
30  bool verbose;
31  std::string settingsFile;
32  std::string imageType;
33 };
34 
35 struct BandpassSettings
36 {
37  float highPass;
38  float lowPass;
39  float alphaLow;
40  float alphaHigh;
41  float speedOfSound;
42 };
43 
44 struct CropSettings
45 {
46  int above;
47  int below;
48  int right;
49  int left;
50  int zStart;
51  int zEnd;
52 };
53 
54 struct ResampleSettings
55 {
56  double spacing;
57  unsigned int dimX;
58 };
59 
60 struct BModeSettings
61 {
63  bool UseLogFilter;
64 };
65 
66 struct ProcessSettings
67 {
68  bool DoBeamforming;
69  bool DoBandpass;
70  bool DoCropping;
71  bool DoResampling;
72  bool DoBmode;
73 };
74 
75 InputParameters parseInput(int argc, char* argv[])
76 {
77  mitkCommandLineParser parser;
78  parser.setCategory("MITK-Photoacoustics");
79  parser.setTitle("Mitk Photoacoustics Beamforming Tool");
80  parser.setDescription("Reads a nrrd file as an input and applies a beamforming method as set with the parameters defined in an additionally provided xml file.");
81  parser.setContributor("Computer Assisted Medical Interventions, DKFZ");
82 
83  parser.setArgumentPrefix("--", "-");
84 
85  parser.beginGroup("Required parameters");
86  parser.addArgument(
87  "inputImage", "i", mitkCommandLineParser::Image,
88  "Input image (mitk::Image)", "input image (.nrrd file)",
89  us::Any(), false, false, false, mitkCommandLineParser::Input);
90  parser.addArgument(
91  "output", "o", mitkCommandLineParser::File,
92  "Output filename", "output image (.nrrd file)",
93  us::Any(), false, false, false, mitkCommandLineParser::Output);
94  parser.addArgument(
95  "settings", "s", mitkCommandLineParser::String,
96  "settings file", "file containing beamforming and other specifications(.xml file)",
97  us::Any(), false);
98  parser.addArgument(
99  "type", "t", mitkCommandLineParser::String,
100  "image type", "Specifies whether to use the PA or US subsection of the xml file must be 'PA' or 'US'",
101  us::Any(), false);
102  parser.endGroup();
103 
104  parser.beginGroup("Optional parameters");
105  parser.addArgument(
106  "verbose", "v", mitkCommandLineParser::Bool,
107  "Verbose Output", "Whether to produce verbose, or rather debug output. (default: false)");
108  parser.endGroup();
109 
110  InputParameters input;
111 
112  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
113  if (parsedArgs.size() == 0)
114  exit(-1);
115 
116  input.verbose = (bool)parsedArgs.count("verbose");
117  MITK_INFO(input.verbose) << "### VERBOSE OUTPUT ENABLED ###";
118 
119  if (parsedArgs.count("inputImage"))
120  {
121  MITK_INFO(input.verbose) << "Reading input image...";
122  input.inputImage = mitk::IOUtil::Load<mitk::Image>(us::any_cast<std::string>(parsedArgs["inputImage"]));
123  MITK_INFO(input.verbose) << "Reading input image...[Done]";
124  }
125  else
126  mitkThrow() << "No input image given.";
127 
128  if (parsedArgs.count("output"))
129  input.outputFilename = us::any_cast<std::string>(parsedArgs["output"]);
130  else
131  mitkThrow() << "No output image path given..";
132 
133  if (parsedArgs.count("settings"))
134  input.settingsFile = us::any_cast<std::string>(parsedArgs["settings"]);
135  else
136  mitkThrow() << "No settings image path given..";
137 
138  if (parsedArgs.count("type"))
139  input.imageType = us::any_cast<std::string>(parsedArgs["type"]);
140  else
141  mitkThrow() << "No settings image type given..";
142 
143  return input;
144 }
145 
146 TiXmlElement* getRootChild(TiXmlElement* root, std::string childName)
147 {
148  for (TiXmlElement* elem = root->FirstChildElement(); elem != NULL; elem = elem->NextSiblingElement())
149  {
150  if (elem->Value() == childName)
151  return elem;
152  }
153  return nullptr;
154 }
155 
156 void ParseXML(std::string xmlFile, InputParameters input, mitk::BeamformingSettings::Pointer *bfSet, BandpassSettings& bandpassSet, CropSettings& cropSet, ResampleSettings& resSet, BModeSettings& bmodeSet, ProcessSettings& processSet)
157 {
158  MITK_INFO << "Loading configuration File \"" << xmlFile << "\"";
159  TiXmlDocument doc(xmlFile);
160  if (!doc.LoadFile())
161  mitkThrow() << "Failed to load settings file \"" << xmlFile << "\" Error: " << doc.ErrorDesc();
162 
163  TiXmlElement* root = doc.FirstChildElement();
164  if (root == NULL)
165  {
166  mitkThrow() << "Failed to load file: No root element.";
167  doc.Clear();
168  }
169 
170  TiXmlElement* elemtop = getRootChild(root, input.imageType);
171  if(elemtop == nullptr)
172  mitkThrow() << "The xml file is wrongly formatted, no element \"" << input.imageType << "\" found";
173 
174  for (TiXmlElement* elem = elemtop->FirstChildElement(); elem != NULL; elem = elem->NextSiblingElement())
175  {
176  std::string elemName = elem->Value();
177  if (elemName == "Beamforming")
178  {
179  processSet.DoBeamforming = std::stoi(elem->Attribute("do"));
180  if (!processSet.DoBeamforming)
181  continue;
182 
183  float SpeedOfSound = std::stof(elem->Attribute("speedOfSoundMeterPerSecond"));
184  float PitchMeter = std::stof(elem->Attribute("pitchMilliMeter"))/1000;
185  float TimeSpacingMicroSecond = (float)(input.inputImage->GetGeometry()->GetSpacing()[1] / 1000000);
186 
187  if (std::stof(elem->Attribute("timeSpacingMicroSecond")) > 0) {
188  TimeSpacingMicroSecond = std::stof(elem->Attribute("timeSpacingMicroSecond"));
189  MITK_INFO << "TIME: " << TimeSpacingMicroSecond;
190  }
191 
192  float Angle = std::stof(elem->Attribute("apodizationAngle"));
193  bool IsPhotoacousticImage = true;
194  if (std::string(elem->Attribute("imageType")).compare("US") == 0) {
195  IsPhotoacousticImage = false;
196  MITK_INFO << "US IMAGE";
197  }
198 
199  unsigned int SamplesPerLine = std::stoi(elem->Attribute("reconstructedYDimension"));
200  unsigned int ReconstructionLines = std::stoi(elem->Attribute("reconstructedXDimension"));
201 
202 
203  float ReconstructionDepth = std::stof(elem->Attribute("reconstructionDepthMeter"));
204  bool UseGPU = std::stoi(elem->Attribute("useGPU"));
205  unsigned int GPUBatchSize = std::stoi(elem->Attribute("GPUBatchSize"));
206 
207  std::string apodizationStr = elem->Attribute("apodizationFunction");
208  mitk::BeamformingSettings::Apodization Apodization = mitk::BeamformingSettings::Apodization::Box;
209  if (apodizationStr == "Box")
210  Apodization = mitk::BeamformingSettings::Apodization::Box;
211  else if (apodizationStr == "Hann")
212  Apodization = mitk::BeamformingSettings::Apodization::Hann;
213  else if (apodizationStr == "Hamm")
214  Apodization = mitk::BeamformingSettings::Apodization::Hamm;
215  else
216  mitkThrow() << "Apodization incorrectly defined in settings";
217 
218  std::string geomStr = elem->Attribute("probeGeometry");
220  if (geomStr == "linear")
222  else if(geomStr == "concave")
223  ProbeGeometry = mitk::BeamformingSettings::ProbeGeometry::Concave;
224  else
225  mitkThrow() << "geometry incorrectly defined in settings";
226 
227  float radius = std::stof(elem->Attribute("radiusMilliMeter"))/1000;
228 
229  unsigned int ApodizationArraySize = ReconstructionLines;
230 
231  std::string algorithmStr = elem->Attribute("algorithm");
232  mitk::BeamformingSettings::BeamformingAlgorithm Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS;
233  if (algorithmStr == "DAS")
234  Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS;
235  else if (algorithmStr == "DMAS")
236  Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS;
237  else if (algorithmStr == "sDMAS")
238  Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS;
239  else
240  mitkThrow() << "Beamforming algorithm incorrectly defined in settings";
241 
243  PitchMeter,
244  SpeedOfSound,
245  TimeSpacingMicroSecond,
246  Angle,
247  IsPhotoacousticImage,
248  SamplesPerLine,
249  ReconstructionLines,
250  input.inputImage->GetDimensions(),
251  ReconstructionDepth,
252  UseGPU,
253  GPUBatchSize,
254  Apodization,
255  ApodizationArraySize,
256  Algorithm,
257  ProbeGeometry,
258  radius
259  );
260  }
261  if (elemName == "Bandpass")
262  {
263  processSet.DoBandpass = std::stoi(elem->Attribute("do"));
264  if (!processSet.DoBandpass)
265  continue;
266 
267  bandpassSet.highPass = std::stof(elem->Attribute("highPass"));
268  bandpassSet.lowPass = std::stof(elem->Attribute("lowPass"));
269  bandpassSet.alphaLow = std::stof(elem->Attribute("alphaLow"));
270  bandpassSet.alphaHigh = std::stof(elem->Attribute("alphaHigh"));
271  }
272  if (elemName == "Cropping")
273  {
274  processSet.DoCropping = std::stoi(elem->Attribute("do"));
275  if (!processSet.DoCropping)
276  continue;
277 
278  cropSet.above = std::stoi(elem->Attribute("cutAbove"));
279  cropSet.below = std::stoi(elem->Attribute("cutBelow"));
280  cropSet.right = std::stoi(elem->Attribute("cutRight"));
281  cropSet.left = std::stoi(elem->Attribute("cutLeft"));
282  cropSet.zStart = std::stoi(elem->Attribute("firstSlice"));
283  cropSet.zEnd = std::stoi(elem->Attribute("cutSlices"));
284  }
285  if (elemName == "Resampling")
286  {
287  processSet.DoResampling = std::stoi(elem->Attribute("do"));
288  if (!processSet.DoResampling)
289  continue;
290 
291  resSet.spacing = std::stod(elem->Attribute("spacing"));
292  resSet.dimX = std::stoi(elem->Attribute("dimX"));
293  }
294  if (elemName == "BMode")
295  {
296  processSet.DoBmode = std::stoi(elem->Attribute("do"));
297  if (!processSet.DoBmode)
298  continue;
299 
300  std::string methodStr = elem->Attribute("method");
301  if (methodStr == "EnvelopeDetection")
302  bmodeSet.method = mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection;
303  else if (methodStr == "Abs")
304  bmodeSet.method = mitk::PhotoacousticFilterService::BModeMethod::Abs;
305  else
306  mitkThrow() << "BMode method incorrectly set in configuration file";
307  bmodeSet.UseLogFilter = (bool)std::stoi(elem->Attribute("useLogFilter"));
308  }
309  }
310 }
311 
312 int main(int argc, char * argv[])
313 {
314  auto input = parseInput(argc, argv);
315 
316  mitk::BeamformingSettings::Pointer bfSettings;
317  BandpassSettings bandpassSettings{5,10,1,1,1540};
318  BModeSettings bmodeSettings{ mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, false };
319  CropSettings cropSettings{ 0,0,0,0,0,0 };
320  ResampleSettings resSettings{ 0.15 , 256};
321  ProcessSettings processSettings{ false, false, false, false, false };
322 
323  MITK_INFO << "Parsing settings XML...";
324  try
325  {
326  ParseXML(input.settingsFile, input, &bfSettings, bandpassSettings, cropSettings, resSettings, bmodeSettings, processSettings);
327  }
328  catch (mitk::Exception e)
329  {
330  MITK_INFO << e;
331  return -1;
332  }
333 
334  MITK_INFO << "Parsing settings XML...[Done]";
335 
336  MITK_INFO(input.verbose) << "Beamforming input image...";
337  mitk::Image::Pointer inputImage = input.inputImage;
338  if (!(inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)"))
339  {
340  // we need a float image, so cast it here
341  MITK_INFO(input.verbose) << "Casting input image to float...";
343  castFilter->SetInput(inputImage);
344  castFilter->Update();
345  inputImage = castFilter->GetOutput();
346  MITK_INFO << inputImage->GetPixelType().GetPixelTypeAsString();
347  MITK_INFO(input.verbose) << "Casting input image to float...[Done]";
348  }
349 
350  mitk::PhotoacousticFilterService::Pointer m_FilterService = mitk::PhotoacousticFilterService::New();
351 
352  mitk::Image::Pointer output = inputImage;
353  if (processSettings.DoBandpass)
354  {
355  MITK_INFO(input.verbose) << "Bandpassing input image...";
356  output = m_FilterService->ApplyBandpassFilter(output, bandpassSettings.highPass*1e6, bandpassSettings.lowPass*1e6, bandpassSettings.alphaHigh, bandpassSettings.alphaLow, 1, bandpassSettings.speedOfSound, true);
357  MITK_INFO(input.verbose) << "Bandpassing input image...[Done]";
358  }
359  if (processSettings.DoBeamforming)
360  {
361  MITK_INFO(input.verbose) << "Beamforming input image...";
362  output = m_FilterService->ApplyBeamforming(output, bfSettings);
363  MITK_INFO(input.verbose) << "Beamforming input image...[Done]";
364  }
365  if (processSettings.DoCropping)
366  {
367  int err;
368  MITK_INFO(input.verbose) << "Applying Crop filter to image...";
369  output = m_FilterService->ApplyCropping(output,
370  cropSettings.above, cropSettings.below, cropSettings.right, cropSettings.left, cropSettings.zStart, cropSettings.zEnd, &err);
371  MITK_INFO(input.verbose) << "Applying Crop filter to image...[Done]";
372  }
373  if (processSettings.DoBmode)
374  {
375  MITK_INFO(input.verbose) << "Applying BModeFilter to image...";
376  output = m_FilterService->ApplyBmodeFilter(output, bmodeSettings.method, bmodeSettings.UseLogFilter);
377  MITK_INFO(input.verbose) << "Applying BModeFilter to image...[Done]";
378  }
379  if (processSettings.DoResampling)
380  {
381  double spacing[3] = { output->GetGeometry()->GetSpacing()[0] * output->GetDimension(0) / resSettings.dimX, resSettings.spacing, output->GetGeometry()->GetSpacing()[2] };
382  MITK_INFO(input.verbose) << "Applying Resample filter to image...";
383  output = m_FilterService->ApplyResampling(output, spacing);
384  if (output->GetDimension(0) != resSettings.dimX)
385  {
386  double dim[3] = { (double)resSettings.dimX, (double)output->GetDimension(1), (double)output->GetDimension(2) };
387  output = m_FilterService->ApplyResamplingToDim(output, dim);
388  }
389  MITK_INFO(input.verbose) << "Applying Resample filter to image...[Done]";
390  }
391 
392  MITK_INFO(input.verbose) << "Saving image...";
393  mitk::IOUtil::Save(output, input.outputFilename);
394  MITK_INFO(input.verbose) << "Saving image...[Done]";
395 
396  MITK_INFO(input.verbose) << "Beamforming input image...[Done]";
397 }
#define MITK_INFO
Definition: mitkLogMacros.h:18
InputParameters parseInput(int argc, char *argv[])
void setContributor(std::string contributor)
void ParseXML(std::string xmlFile, InputParameters input, mitk::BeamformingSettings::Pointer *bfSet, BandpassSettings &bandpassSet, CropSettings &cropSet, ResampleSettings &resSet, BModeSettings &bmodeSet, ProcessSettings &processSet)
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
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)
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
std::string outputFilename
Definition: MitkMCxyz.cpp:626
ProbeGeometry
Available geometries for Probes:
An object of this class represents an exception of MITK. Please don&#39;t instantiate exceptions manually...
Definition: mitkException.h:45
bool compare(std::pair< double, int > i, std::pair< double, int > j)
int main(int argc, char *argv[])
#define mitkThrow()
bool verbose(false)
Definition: usAny.h:163
BeamformingAlgorithm
Available beamforming algorithms:
static Pointer New(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int *inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, Apodization apod, unsigned int apodizationArraySize, BeamformingAlgorithm algorithm, ProbeGeometry geometry, float probeRadius)
void setCategory(std::string category)
Apodization
Available apodization functions:
void setArgumentPrefix(const std::string &longPrefix, const std::string &shortPrefix)
static void Save(const mitk::BaseData *data, const std::string &path, bool setPathProperty=false)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:774
void setTitle(std::string title)
void setDescription(std::string description)
TiXmlElement * getRootChild(TiXmlElement *root, std::string childName)
BModeMethod
Defines the methods for the B-Mode filter Currently implemented are an Envelope Detection filter and ...
void beginGroup(const std::string &description)