Medical Imaging Interaction Toolkit  2018.4.99-1640525a
Medical Imaging Interaction Toolkit
mitkPhotoacousticOCLBeamformingFilter.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 #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN
14 
16 #include "usServiceReference.h"
17 
19  m_PixelCalculation(NULL),
20  m_inputSlices(1),
21  m_Conf(settings),
22  m_InputImage(mitk::Image::New()),
23  m_ApodizationBuffer(nullptr),
24  m_DelaysBuffer(nullptr),
25  m_UsedLinesBuffer(nullptr),
26  m_ElementHeightsBuffer(nullptr),
27  m_ElementPositionsBuffer(nullptr)
28 {
29  MITK_INFO << "Instantiating OCL beamforming Filter...";
30  this->AddSourceFile("DAS.cl");
31  this->AddSourceFile("DMAS.cl");
32  this->AddSourceFile("sDMAS.cl");
33  this->m_FilterID = "OpenCLBeamformingFilter";
34 
35  this->Initialize();
36 
37  unsigned int dim[] = { 128, 2048, 2 };
38 
39  m_InputImage->Initialize(mitk::MakeScalarPixelType<float>(), 3, dim);
40 
41  m_ChunkSize[0] = 128;
42  m_ChunkSize[1] = 128;
43  m_ChunkSize[2] = 8;
44 
45  m_UsedLinesCalculation = mitk::OCLUsedLinesCalculation::New(m_Conf);
46  m_DelayCalculation = mitk::OCLDelayCalculation::New(m_Conf);
47  MITK_INFO << "Instantiating OCL beamforming Filter...[Done]";
48 }
49 
51 {
52  if (this->m_PixelCalculation)
53  {
54  clReleaseKernel(m_PixelCalculation);
55  }
56 
57  if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer);
58  if (m_ElementHeightsBuffer) clReleaseMemObject(m_ElementHeightsBuffer);
59  if (m_ElementPositionsBuffer) clReleaseMemObject(m_ElementPositionsBuffer);
60 }
61 
63 {
64  //Check if context & program available
65  if (!this->Initialize())
66  {
69 
70  // clean-up also the resources
71  resources->InvalidateStorage();
72  mitkThrow() << "Filter is not initialized. Cannot update.";
73  }
74  else {
75  // Execute
76  this->Execute();
77  }
78 }
79 
80 void mitk::PhotoacousticOCLBeamformingFilter::UpdateDataBuffers()
81 {
84  cl_context gpuContext = resources->GetContext();
85 
86  //Initialize the Output
87  try
88  {
89  size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() *
90  (size_t)m_inputSlices;
91  m_OutputDim[0] = m_Conf->GetReconstructionLines();
92  m_OutputDim[1] = m_Conf->GetSamplesPerLine();
93  m_OutputDim[2] = m_inputSlices;
94  this->InitExec(this->m_PixelCalculation, m_OutputDim, outputSize, sizeof(float));
95  }
96  catch (const mitk::Exception& e)
97  {
98  MITK_ERROR << "Caught exception while initializing filter: " << e.what();
99  return;
100  }
101 
102  cl_int clErr = 0;
103  MITK_DEBUG << "Updating GPU Buffers for new configuration";
104 
105  // create the apodisation buffer
106  if (m_Apodisation == nullptr)
107  {
108  MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation.";
109  m_Apodisation = new float[1]{ 1 };
110  m_ApodArraySize = 1;
111  }
112 
113  if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer);
114 
115  this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, const_cast<float*>(m_Apodisation), &clErr);
116  CHECK_OCL_ERR(clErr);
117 
118  if (m_ElementHeightsBuffer) clReleaseMemObject(m_ElementHeightsBuffer);
119  this->m_ElementHeightsBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_Conf->GetTransducerElements(), const_cast<float*>(m_Conf->GetElementHeights()), &clErr);
120  CHECK_OCL_ERR(clErr);
121 
122  if (m_ElementPositionsBuffer) clReleaseMemObject(m_ElementPositionsBuffer);
123  this->m_ElementPositionsBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_Conf->GetTransducerElements(), const_cast<float*>(m_Conf->GetElementPositions()), &clErr);
124  CHECK_OCL_ERR(clErr);
125  // calculate used lines
126 
127  m_UsedLinesCalculation->SetElementPositionsBuffer(m_ElementPositionsBuffer);
128  m_UsedLinesCalculation->SetElementHeightsBuffer(m_ElementHeightsBuffer);
129  m_UsedLinesCalculation->Update();
130  m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer();
131 
132  // calculate the Delays
133  m_DelayCalculation->SetInputs(m_UsedLinesBuffer);
134  m_DelayCalculation->Update();
135 
136  m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer();
137 }
138 
139 void mitk::PhotoacousticOCLBeamformingFilter::Execute()
140 {
141  cl_int clErr = 0;
142  UpdateDataBuffers();
143  if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear)
144  {
145  unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines();
146  unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine();
147 
148  clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer));
149  clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer));
150  clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer));
151  clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize));
152  clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0]));
153  clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1]));
154  clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(m_inputSlices));
155  clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(reconstructionLines));
156  clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine));
157  }
158  else
159  {
160  int reconstructionLines = this->m_Conf->GetReconstructionLines();
161  int samplesPerLine = this->m_Conf->GetSamplesPerLine();
162  float totalSamples_i = (float)(m_Conf->GetReconstructionDepth()) / (float)(m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing());
163  totalSamples_i = totalSamples_i <= m_Conf->GetInputDim()[1] ? totalSamples_i : m_Conf->GetInputDim()[1];
164  float horizontalExtent = m_Conf->GetHorizontalExtent();
165  float mult = 1 / (this->m_Conf->GetTimeSpacing() * this->m_Conf->GetSpeedOfSound());
166  char isPAImage = (char)m_Conf->GetIsPhotoacousticImage();
167 
168  clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_ElementHeightsBuffer));
169  clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_ElementPositionsBuffer));
170  clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer));
171  clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize));
172  clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0]));
173  clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1]));
174  clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_int), &(m_inputSlices));
175  clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_int), &(reconstructionLines));
176  clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_int), &(samplesPerLine));
177  clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_float), &(totalSamples_i));
178  clErr |= clSetKernelArg(this->m_PixelCalculation, 12, sizeof(cl_float), &(horizontalExtent));
179  clErr |= clSetKernelArg(this->m_PixelCalculation, 13, sizeof(cl_float), &(mult));
180  clErr |= clSetKernelArg(this->m_PixelCalculation, 14, sizeof(cl_char), &(isPAImage));
181  clErr |= clSetKernelArg(this->m_PixelCalculation, 15, sizeof(cl_mem), &(this->m_UsedLinesBuffer));
182  }
183  // execute the filter on a 2D/3D NDRange
184  if (m_OutputDim[2] == 1 || m_ChunkSize[2] == 1)
185  {
186  if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, m_inputSlices, 50))
187  mitkThrow() << "openCL Error when executing Kernel";
188  }
189  else
190  {
191  if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 3, m_ChunkSize, m_inputSlices, 50))
192  mitkThrow() << "openCL Error when executing Kernel";
193  }
194 
195  // signalize the GPU-side data changed
196  m_Output->Modified(GPU_DATA);
197 }
198 
199 us::Module *mitk::PhotoacousticOCLBeamformingFilter::GetModule()
200 {
201  return us::GetModuleContext()->GetModule();
202 }
203 
204 bool mitk::PhotoacousticOCLBeamformingFilter::Initialize()
205 {
206  bool buildErr = true;
207  cl_int clErr = 0;
208 
209  if (OclFilter::Initialize())
210  {
211  if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear)
212  {
213  switch (m_Conf->GetAlgorithm())
214  {
215  case BeamformingSettings::BeamformingAlgorithm::DAS:
216  {
217  MITK_INFO << "DAS bf";
218  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr);
219  break;
220  }
221  case BeamformingSettings::BeamformingAlgorithm::DMAS:
222  {
223  MITK_INFO << "DMAS bf";
224  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr);
225  break;
226  }
227  case BeamformingSettings::BeamformingAlgorithm::sDMAS:
228  {
229  MITK_INFO << "sDMAS bf";
230  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS", &clErr);
231  break;
232  }
233  default:
234  {
235  MITK_INFO << "No beamforming algorithm specified, setting to DAS";
236  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr);
237  break;
238  }
239  }
240  }
241  else
242  {
243  switch (m_Conf->GetAlgorithm())
244  {
245  case BeamformingSettings::BeamformingAlgorithm::DAS:
246  {
247  MITK_INFO << "DAS bf";
248  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS_g", &clErr);
249  break;
250  }
251  case BeamformingSettings::BeamformingAlgorithm::DMAS:
252  {
253  MITK_INFO << "DMAS bf";
254  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS_g", &clErr);
255  break;
256  }
257  case BeamformingSettings::BeamformingAlgorithm::sDMAS:
258  {
259  MITK_INFO << "sDMAS bf";
260  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS_g", &clErr);
261  break;
262  }
263  default:
264  {
265  MITK_INFO << "No beamforming algorithm specified, setting to DAS";
266  this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS_g", &clErr);
267  break;
268  }
269  }
270  }
271  buildErr |= CHECK_OCL_ERR(clErr);
272  }
273 
274  CHECK_OCL_ERR(clErr);
275 
276  return (OclFilter::IsInitialized() && buildErr);
277 }
278 
279 void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image)
280 {
282  m_InputImage = image;
283  m_inputSlices = image->GetDimension(2);
284 }
285 
286 void mitk::PhotoacousticOCLBeamformingFilter::SetInput(void* data, unsigned int* dimensions, unsigned int BpE)
287 {
288  OclDataSetToDataSetFilter::SetInput(data, dimensions[0] * dimensions[1] * dimensions[2], BpE);
289 }
290 
291 mitk::Image::Pointer mitk::PhotoacousticOCLBeamformingFilter::GetOutputAsImage()
292 {
293  mitk::Image::Pointer outputImage = mitk::Image::New();
294 
295  if (m_Output->IsModified(GPU_DATA))
296  {
297  void* pData = m_Output->TransferDataToCPU(m_CommandQue);
298 
299  const unsigned int dimension = 3;
300  unsigned int dimensions[3] = { (unsigned int)m_OutputDim[0], (unsigned int)m_OutputDim[1], (unsigned int)m_OutputDim[2] };
301 
302  const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry();
303 
304  MITK_DEBUG << "Creating new MITK Image.";
305 
306  outputImage->Initialize(this->GetOutputType(), dimension, dimensions);
307  outputImage->SetSpacing(p_slg->GetSpacing());
308  outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory);
309  free(pData);
310  }
311 
312  MITK_DEBUG << "Image Initialized.";
313 
314  return outputImage;
315 }
316 
317 void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput()
318 {
320 }
321 #endif
#define CHECK_OCL_ERR(_er)
Definition: mitkOclUtils.h:21
virtual void InvalidateStorage()=0
Remove all invalid (=do not compile) programs from the internal storage.
ServiceReferenceU GetServiceReference(const std::string &clazz)
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define GPU_DATA
#define MITK_ERROR
Definition: mitkLogMacros.h:20
static void Update(vtkPolyData *)
Definition: mitkSurface.cpp:31
bool Initialize()
Initialize all necessary parts of the filter.
#define MITK_DEBUG
Definition: mitkLogMacros.h:22
DataCollection - Class to facilitate loading/accessing structured data.
virtual bool IsInitialized()
Returns true if the initialization was successfull.
void * GetService(const ServiceReferenceBase &reference)
class ITK_EXPORT Image
void SetInput(mitk::OclDataSet::Pointer DataSet)
SetInput SetInput Set the input DataSet (as mitk::OclDataSet).
void * GetOutput()
Returns an pointer to the filtered data.
An object of this class represents an exception of MITK. Please don&#39;t instantiate exceptions manually...
Definition: mitkException.h:45
#define mitkThrow()
Module * GetModule() const
mitk::Image::Pointer image
static Pointer New()
virtual cl_context GetContext() const =0
Returns a valid OpenCL Context (if applicable) or nullptr if none present.
static ModuleContext * GetModuleContext()
Returns the module context of the calling module.
Declaration of the OpenCL Resources micro-service.