Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkBeamformingUtils.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 "mitkBeamformingUtils.h"
22 
24 {
25 }
26 
28 {
29 }
30 
32 {
33  float* ApodWindow = new float[samples];
34 
35  for (int n = 0; n < samples; ++n)
36  {
37  ApodWindow[n] = (1 - cos(2 * itk::Math::pi * n / (samples - 1))) / 2;
38  }
39 
40  return ApodWindow;
41 }
42 
44 {
45  float* ApodWindow = new float[samples];
46 
47  for (int n = 0; n < samples; ++n)
48  {
49  ApodWindow[n] = 0.54 - 0.46*cos(2 * itk::Math::pi*n / (samples - 1));
50  }
51 
52  return ApodWindow;
53 }
54 
56 {
57  float* ApodWindow = new float[samples];
58 
59  for (int n = 0; n < samples; ++n)
60  {
61  ApodWindow[n] = 1;
62  }
63 
64  return ApodWindow;
65 }
66 
67 unsigned short* mitk::BeamformingUtils::MinMaxLines(const mitk::BeamformingSettings::Pointer config)
68 {
69  unsigned int outputL = (unsigned int)config->GetReconstructionLines();
70  unsigned int outputS = (unsigned int)config->GetSamplesPerLine();
71 
72  unsigned short* dDest = new unsigned short[outputL * outputS * 2];
73 
74  unsigned int inputL = (unsigned int)config->GetInputDim()[0];
75 
76  float horizontalExtent = config->GetHorizontalExtent();
77  float verticalExtent = config->GetReconstructionDepth();
78 
79  float partMult = (tan(config->GetAngle() / 360 * 2 * itk::Math::pi) *
80  ((config->GetSpeedOfSound() * config->GetTimeSpacing())) /
81  (config->GetPitchInMeters() * config->GetTransducerElements())) * inputL;
82  float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing());
83 
84  totalSamples_i = totalSamples_i <= config->GetInputDim()[1] ? totalSamples_i : config->GetInputDim()[1];
85 
86  if ((int)config->GetGeometry() == 0) // if this is raw data from a linear probe geometry
87  {
88  for (unsigned int x = 0; x < outputL; ++x)
89  {
90  for (unsigned int y = 0; y < outputS; ++y)
91  {
92  float l_i = (float)x / outputL * inputL;
93  float s_i = (float)y / (float)outputS * totalSamples_i;
94 
95  float part = partMult * s_i;
96  if (part < 1)
97  part = 1;
98  unsigned short maxLine = std::min((l_i + part) + 1, (float)inputL);
99  unsigned short minLine = std::max((l_i - part), 0.0f);
100 
101  dDest[y * 2 * outputL + 2 * x] = (unsigned short)minLine; //minLine
102  dDest[y * 2 * outputL + 2 * x + 1] = (unsigned short)maxLine; //maxLine
103  }
104  }
105 
106  }
107  else // if this is *not* raw data from a linear probe geometry (currently meaning its a concave geometry)
108  {
109  float probeRadius = config->GetProbeRadius();
110  float* elementHeights = config->GetElementHeights();
111  float* elementPositions = config->GetElementPositions();
112 
113  float sin_deg = std::sin(config->GetAngle() / 360 * 2 * itk::Math::pi);
114 
115  float x_center_pos = horizontalExtent / 2.0;
116  float y_center_pos = probeRadius;
117 
118  for (unsigned int x = 0; x < outputL; ++x)
119  {
120  for (unsigned int y = 0; y < outputS; ++y)
121  {
122  float x_cm = (float)x / outputL * horizontalExtent; // x*Xspacing
123  float y_cm = (float)y / (float)outputS * verticalExtent; // y*Yspacing
124 
125  unsigned int maxLine = inputL;
126  unsigned int minLine = 0;
127 
128  for (unsigned int l_s = minLine; l_s <= inputL; l_s += 1)
129  {
130  float x_sensor_pos = elementPositions[l_s];
131  float y_sensor_pos = elementHeights[l_s];
132 
133  float distance_sensor_target = sqrt((x_cm - x_sensor_pos)*(x_cm - x_sensor_pos)
134  + (y_cm - y_sensor_pos)*(y_cm - y_sensor_pos));
135 
136  // solving line equation
137  float center_to_sensor_a = y_sensor_pos - y_center_pos;
138  float center_to_sensor_b = x_center_pos - x_sensor_pos;
139  float center_to_sensor_c = -(center_to_sensor_a * x_center_pos + center_to_sensor_b * y_center_pos);
140  float distance_to_sensor_direction = std::fabs((center_to_sensor_a * x_cm
141  + center_to_sensor_b * y_cm
142  + center_to_sensor_c)) /
143  (sqrt(center_to_sensor_a*center_to_sensor_a + center_to_sensor_b*center_to_sensor_b));
144 
145  if (distance_to_sensor_direction < sin_deg*distance_sensor_target)
146  {
147  minLine = l_s;
148  break;
149  }
150  }
151 
152  for (unsigned int l_s = maxLine - 1; l_s > minLine; l_s -= 1) // start with maxline-1 otherwise elementPositions[] will go out of range
153  {
154  float x_sensor_pos = elementPositions[l_s];
155  float y_sensor_pos = elementHeights[l_s];
156 
157  float distance_sensor_target = sqrt((x_cm - x_sensor_pos)*(x_cm - x_sensor_pos)
158  + (y_cm - y_sensor_pos)*(y_cm - y_sensor_pos));
159 
160  // solving line equation
161  float center_to_sensor_a = y_sensor_pos - y_center_pos;
162  float center_to_sensor_b = x_center_pos - x_sensor_pos;
163  float center_to_sensor_c = -(center_to_sensor_a * x_center_pos + center_to_sensor_b * y_center_pos);
164  float distance_to_sensor_direction = std::fabs((center_to_sensor_a * x_cm
165  + center_to_sensor_b * y_cm
166  + center_to_sensor_c)) /
167  (sqrt(center_to_sensor_a*center_to_sensor_a + center_to_sensor_b*center_to_sensor_b));
168 
169  if (distance_to_sensor_direction < sin_deg*distance_sensor_target)
170  {
171  maxLine = l_s;
172  break;
173  }
174  }
175  dDest[y * 2 * outputL + 2 * x] = (unsigned short)minLine; //minLine
176  dDest[y * 2 * outputL + 2 * x + 1] = (unsigned short)maxLine; //maxLine
177  }
178  }
179  }
180  return dDest;
181 }
182 
184  float* input, float* output, float inputDim[2], float outputDim[2],
185  const short& line, const mitk::BeamformingSettings::Pointer config)
186 {
187  const float* apodisation = config->GetApodizationFunction();
188  const short apodArraySize = config->GetApodizationArraySize();
189 
190  const float* elementHeights = config->GetElementHeights();
191  const float* elementPositions = config->GetElementPositions();
192 
193  float& inputS = inputDim[1];
194  float& inputL = inputDim[0];
195 
196  float& outputS = outputDim[1];
197  float& outputL = outputDim[0];
198 
199  short AddSample = 0;
200  short maxLine = 0;
201  short minLine = 0;
202  float l_p = 0;
203  float s_i = 0;
204 
205  float apod_mult = 1;
206 
207  short usedLines = (maxLine - minLine);
208 
209  float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing());
210  totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS;
211 
212  l_p = (float)line / outputL * config->GetHorizontalExtent();
213 
214  for (short sample = 0; sample < outputS; ++sample)
215  {
216  s_i = (float)sample / outputS * totalSamples_i;
217 
218  minLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line];
219  maxLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line + 1];
220  usedLines = (maxLine - minLine);
221 
222  apod_mult = (float)apodArraySize / (float)usedLines;
223 
224  for (short l_s = minLine; l_s < maxLine; ++l_s)
225  {
226  AddSample = (int)sqrt(
227  pow(s_i-elementHeights[l_s]/(config->GetSpeedOfSound()*config->GetTimeSpacing()), 2)
228  +
229  pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2)
230  ) + (1 - config->GetIsPhotoacousticImage())*s_i;
231  if (AddSample < inputS && AddSample >= 0)
232  output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] *
233  apodisation[(short)((l_s - minLine)*apod_mult)];
234  else
235  --usedLines;
236  }
237  output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines;
238  }
239 }
240 
242  float* input, float* output, float inputDim[2], float outputDim[2],
243  const short& line, const mitk::BeamformingSettings::Pointer config)
244 {
245  const float* apodisation = config->GetApodizationFunction();
246  const short apodArraySize = config->GetApodizationArraySize();
247 
248  const float* elementHeights = config->GetElementHeights();
249  const float* elementPositions = config->GetElementPositions();
250 
251  float& inputS = inputDim[1];
252  float& inputL = inputDim[0];
253 
254  float& outputS = outputDim[1];
255  float& outputL = outputDim[0];
256 
257  short maxLine = 0;
258  short minLine = 0;
259  float l_p = 0;
260  float s_i = 0;
261 
262  float apod_mult = 1;
263 
264  float mult = 0;
265 
266  short usedLines = (maxLine - minLine);
267 
268  float totalSamples_i = (float)(config->GetReconstructionDepth()) /
269  (float)(config->GetSpeedOfSound() * config->GetTimeSpacing());
270  totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS;
271 
272  l_p = (float)line / outputL * config->GetHorizontalExtent();
273 
274  for (short sample = 0; sample < outputS; ++sample)
275  {
276  s_i = (float)sample / outputS * totalSamples_i;
277 
278  minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line];
279  maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1];
280  usedLines = (maxLine - minLine);
281 
282  apod_mult = (float)apodArraySize / (float)usedLines;
283 
284  //calculate the AddSamples beforehand to save some time
285  short* AddSample = new short[maxLine - minLine];
286  for (short l_s = 0; l_s < maxLine - minLine; ++l_s)
287  {
288  AddSample[l_s] = (int)sqrt(
289  pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2)
290  +
291  pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 2)
292  ) + (1 - config->GetIsPhotoacousticImage())*s_i;
293  }
294 
295  float s_1 = 0;
296  float s_2 = 0;
297 
298  for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1)
299  {
300  if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0)
301  {
302  for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2)
303  {
304  if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0)
305  {
306  s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL];
307  s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL];
308 
309  mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)];
310  output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0));
311  }
312  }
313  }
314  else
315  --usedLines;
316  }
317 
318  output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1));
319 
320  delete[] AddSample;
321  }
322 }
323 
325  float* input, float* output, float inputDim[2], float outputDim[2],
326  const short& line, const mitk::BeamformingSettings::Pointer config)
327 {
328  const float* apodisation = config->GetApodizationFunction();
329  const short apodArraySize = config->GetApodizationArraySize();
330 
331  const float* elementHeights = config->GetElementHeights();
332  const float* elementPositions = config->GetElementPositions();
333 
334  float& inputS = inputDim[1];
335  float& inputL = inputDim[0];
336 
337  float& outputS = outputDim[1];
338  float& outputL = outputDim[0];
339 
340  short maxLine = 0;
341  short minLine = 0;
342  float l_p = 0;
343  float s_i = 0;
344 
345  float apod_mult = 1;
346 
347  float mult = 0;
348 
349  short usedLines = (maxLine - minLine);
350 
351  float totalSamples_i = (float)(config->GetReconstructionDepth()) /
352  (float)(config->GetSpeedOfSound() * config->GetTimeSpacing());
353  totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS;
354 
355  l_p = (float)line / outputL * config->GetHorizontalExtent();
356 
357  for (short sample = 0; sample < outputS; ++sample)
358  {
359  s_i = (float)sample / outputS * totalSamples_i;
360 
361  minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line];
362  maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1];
363  usedLines = (maxLine - minLine);
364 
365  apod_mult = (float)apodArraySize / (float)usedLines;
366 
367  //calculate the AddSamples beforehand to save some time
368  short* AddSample = new short[maxLine - minLine];
369  for (short l_s = 0; l_s < maxLine - minLine; ++l_s)
370  {
371  AddSample[l_s] = (int)sqrt(
372  pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2)
373  +
374  pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 2)
375  ) + (1 - config->GetIsPhotoacousticImage())*s_i;
376  }
377 
378  float s_1 = 0;
379  float s_2 = 0;
380  float sign = 0;
381 
382  for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1)
383  {
384  if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0)
385  {
386  s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL];
387  sign += s_1;
388 
389  for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2)
390  {
391  if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0)
392  {
393  s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL];
394 
395  mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)];
396  output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0));
397  }
398  }
399  }
400  else
401  --usedLines;
402  }
403 
404  output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0));
405 
406  delete[] AddSample;
407  }
408 }
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
static float * BoxFunction(int samples)
Function to create a Box apodization window.
static float * VonHannFunction(int samples)
Pointer holding the Von-Hann apodization window for beamforming.
static unsigned short * MinMaxLines(const mitk::BeamformingSettings::Pointer config)
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...
static float * HammFunction(int samples)
Function to create a Hamming apodization window.
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...
static T max(T x, T y)
Definition: svm.cpp:56
static T min(T x, T y)
Definition: svm.cpp:53