Medical Imaging Interaction Toolkit  2018.4.99-6ca56567
Medical Imaging Interaction Toolkit
mitkBeamformingFilter.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 "mitkProperties.h"
14 #include "mitkImageReadAccessor.h"
15 #include <algorithm>
16 #include <itkImageIOBase.h>
17 #include <chrono>
18 #include <thread>
19 #include <itkImageIOBase.h>
20 #include "mitkImageCast.h"
21 #include "mitkBeamformingFilter.h"
22 #include "mitkBeamformingUtils.h"
23 
24 mitk::BeamformingFilter::BeamformingFilter(mitk::BeamformingSettings::Pointer settings) :
25  m_OutputData(nullptr),
26  m_InputData(nullptr),
27  m_Conf(settings)
28 {
29  MITK_INFO << "Instantiating BeamformingFilter...";
30  this->SetNumberOfIndexedInputs(1);
31  this->SetNumberOfRequiredInputs(1);
32 
33  m_ProgressHandle = [](int, std::string) {};
34 #if defined(PHOTOACOUSTICS_USE_GPU)
36 #else
38 #endif
39 
40  MITK_INFO << "Instantiating BeamformingFilter...[Done]";
41 }
42 
43 void mitk::BeamformingFilter::SetProgressHandle(std::function<void(int, std::string)> progressHandle)
44 {
45  m_ProgressHandle = progressHandle;
46 }
47 
49 {
50  MITK_INFO << "Destructed BeamformingFilter";
51 }
52 
54 {
55  Superclass::GenerateInputRequestedRegion();
56 
57  mitk::Image* output = this->GetOutput();
58  mitk::Image* input = const_cast<mitk::Image *> (this->GetInput());
59  if (!output->IsInitialized())
60  {
61  return;
62  }
63 
65 }
66 
68 {
69  mitk::Image::ConstPointer input = this->GetInput();
70  mitk::Image::Pointer output = this->GetOutput();
71 
72  if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
73  return;
74 
75  mitk::Vector3D spacing;
76  spacing[0] = m_Conf->GetHorizontalExtent() / m_Conf->GetReconstructionLines() * 1000;
77  float desiredYSpacing = m_Conf->GetReconstructionDepth() * 1000 / m_Conf->GetSamplesPerLine();
78  float maxYSpacing = m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing() * input->GetDimension(1) / m_Conf->GetSamplesPerLine() * 1000;
79  spacing[1] = desiredYSpacing < maxYSpacing ? desiredYSpacing : maxYSpacing;
80  spacing[2] = 1;
81 
82  unsigned int dim[] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), input->GetDimension(2)};
83  output->Initialize(mitk::MakeScalarPixelType<float>(), 3, dim);
84  output->GetGeometry()->SetSpacing(spacing);
85  output->GetGeometry()->Modified();
86  output->SetPropertyList(input->GetPropertyList()->Clone());
87 
89 }
90 
92 {
93  mitk::Image::Pointer input = this->GetInput();
94  if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)"))
95  {
96  MITK_ERROR << "Pixel type of input needs to be float for this filter to work.";
97  mitkThrow() << "Pixel type of input needs to be float for this filter to work.";
98  }
99 
101  mitk::Image::Pointer output = this->GetOutput();
102 
103  if (!output->IsInitialized())
104  return;
105 
106  auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance...
107 
108  if (!m_Conf->GetUseGPU())
109  {
110  int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1;
111  // the interval at which we update the gui progress bar
112 
113  float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) };
114  float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) };
115 
116  for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied
117  {
118  mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i));
119  m_InputData = (float*)inputReadAccessor.GetData();
120 
121  m_OutputData = new float[m_Conf->GetReconstructionLines()*m_Conf->GetSamplesPerLine()];
122 
123  // fill the image with zeros
124  for (int l = 0; l < outputDim[0]; ++l)
125  {
126  for (int s = 0; s < outputDim[1]; ++s)
127  {
128  m_OutputData[l*(short)outputDim[1] + s] = 0;
129  }
130  }
131 
132  std::thread *threads = new std::thread[(short)outputDim[0]];
133 
134  // every line will be beamformed in a seperate thread
135  if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DAS)
136  {
137  for (short line = 0; line < outputDim[0]; ++line)
138  {
139  threads[line] = std::thread(&BeamformingUtils::DASSphericalLine, m_InputData,
140  m_OutputData, inputDim, outputDim, line, m_Conf);
141  }
142  }
143  else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DMAS)
144  {
145  for (short line = 0; line < outputDim[0]; ++line)
146  {
147  threads[line] = std::thread(&BeamformingUtils::DMASSphericalLine, m_InputData,
148  m_OutputData, inputDim, outputDim, line, m_Conf);
149  }
150  }
151  else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::sDMAS)
152  {
153  for (short line = 0; line < outputDim[0]; ++line)
154  {
155  threads[line] = std::thread(&BeamformingUtils::sDMASSphericalLine, m_InputData,
156  m_OutputData, inputDim, outputDim, line, m_Conf);
157  }
158  }
159  // wait for all lines to finish
160  for (short line = 0; line < outputDim[0]; ++line)
161  {
162  threads[line].join();
163  }
164 
165  output->SetSlice(m_OutputData, i);
166 
167  if (i % progInterval == 0)
168  m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction");
169 
170  delete[] m_OutputData;
171  m_OutputData = nullptr;
172  m_InputData = nullptr;
173  }
174  }
175 #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN
176  else
177  {
178  try
179  {
180  // first, we check whether the data is float, other formats are unsupported
181  if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)"))
182  {
183  MITK_ERROR << "Pixel type is not float, abort";
184  mitkThrow() << "Pixel type is not float, abort";
185  }
186 
187  unsigned long availableMemory = m_BeamformingOclFilter->GetDeviceMemory();
188 
189  unsigned int batchSize = m_Conf->GetGPUBatchSize();
190  unsigned int batches = (unsigned int)((float)input->GetDimension(2) / batchSize) + (input->GetDimension(2) % batchSize > 0);
191 
192  unsigned int batchDim[] = { input->GetDimension(0), input->GetDimension(1), batchSize };
193  unsigned int batchDimLast[] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % batchSize };
194 
195  // the following safeguard is probably only needed for absurdly small GPU memory
196  if((unsigned long)batchSize *
197  ((unsigned long)(batchDim[0] * batchDim[1]) * 4 + // single input image (float)
198  (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 4) // single output image (float)
199  > availableMemory -
200  (unsigned long)(m_Conf->GetReconstructionLines() / 2 * m_Conf->GetSamplesPerLine()) * 2 - // Delays buffer (unsigned short)
201  (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 3 * 2 - // UsedLines buffer (unsigned short)
202  50 * 1024 * 1024)// 50 MB buffer for local data, system purposes etc
203  {
204  MITK_ERROR << "device memory too small for GPU beamforming; try decreasing the batch size";
205  return;
206  }
207 
208  mitk::ImageReadAccessor copy(input);
209 
210  for (unsigned int i = 0; i < batches; ++i)
211  {
212  m_ProgressHandle(100.f * (float)i / (float)batches, "performing reconstruction");
213 
214  mitk::Image::Pointer inputBatch = mitk::Image::New();
215  unsigned int num_Slices = 1;
216  if (i == batches - 1 && (input->GetDimension(2) % batchSize > 0))
217  {
218  inputBatch->Initialize(mitk::MakeScalarPixelType<float>(), 3, batchDimLast);
219  num_Slices = batchDimLast[2];
220  }
221  else
222  {
223  inputBatch->Initialize(mitk::MakeScalarPixelType<float>(), 3, batchDim);
224  num_Slices = batchDim[2];
225  }
226 
227  inputBatch->SetSpacing(input->GetGeometry()->GetSpacing());
228 
229  inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i]));
230 
231  m_BeamformingOclFilter->SetApodisation(m_Conf->GetApodizationFunction(), m_Conf->GetApodizationArraySize());
232  m_BeamformingOclFilter->SetInput(inputBatch);
233  m_BeamformingOclFilter->Update();
234 
235  void* out = m_BeamformingOclFilter->GetOutput();
236 
237  for (unsigned int slice = 0; slice < num_Slices; ++slice)
238  {
239  output->SetImportSlice(
240  &(((float*)out)[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * slice]),
241  batchSize * i + slice, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory);
242  }
243  free(out);
244  }
245  }
246  catch (mitk::Exception &e)
247  {
248  std::string errorMessage = "Caught unexpected exception ";
249  errorMessage.append(e.what());
250  MITK_ERROR << errorMessage;
251 
252  float* dummyData = new float[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * m_Conf->GetInputDim()[2]];
253  output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory);
254  }
255  }
256 #endif
257  m_TimeOfHeaderInitialization.Modified();
258 
259  auto end = std::chrono::high_resolution_clock::now();
260  MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin).count()) / 1000000 << "ms" << std::endl;
261 }
void SetRequestedRegionToLargestPossibleRegion() override
static void sDMASSphericalLine(float *input, float *output, float inputDim[2], float outputDim[2], const short &line, const mitk::BeamformingSettings::Pointer config)
Function to perform beamforming on CPU for a single line, using signed DMAS and spherical delay...
static char * line
Definition: svm.cpp:2870
std::function< void(int, std::string)> m_ProgressHandle
The std::function<void(int, std::string)>, through which progress of the currently updating filter is...
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define MITK_ERROR
Definition: mitkLogMacros.h:20
itk::TimeStamp m_TimeOfHeaderInitialization
static void DASSphericalLine(float *input, float *output, float inputDim[2], float outputDim[2], const short &line, const mitk::BeamformingSettings::Pointer config)
Function to perform beamforming on CPU for a single line, using DAS and spherical delay...
void SetProgressHandle(std::function< void(int, std::string)> progressHandle)
Sets a callback for progress checking.
An object of this class represents an exception of MITK. Please don&#39;t instantiate exceptions manually...
Definition: mitkException.h:45
void GenerateInputRequestedRegion() override
#define mitkThrow()
Image class for storing images.
Definition: mitkImage.h:72
static void DMASSphericalLine(float *input, float *output, float inputDim[2], float outputDim[2], const short &line, const mitk::BeamformingSettings::Pointer config)
Function to perform beamforming on CPU for a single line, using DMAS and spherical delay...
BeamformingSettings::Pointer m_Conf
Current configuration set.
static Pointer New()
InputImageType * GetInput(void)
virtual bool IsInitialized() const
Check whether the data has been initialized, i.e., at least the Geometry and other header data has be...
mitk::PhotoacousticOCLBeamformingFilter::Pointer m_BeamformingOclFilter
Pointer to the GPU beamforming filter class; for performance reasons the filter is initialized within...
OutputType * GetOutput()
Get the output data of this image source object.
ImageReadAccessor class to get locked read access for a particular image part.
BeamformingFilter(mitk::BeamformingSettings::Pointer settings)
const void * GetData() const
Gives const access to the data.
void GenerateOutputInformation() override