Medical Imaging Interaction Toolkit  2018.4.99-a3d2e8fb
Medical Imaging Interaction Toolkit
mitkUSDiPhASImageSource.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 // std dependencies
14 #include <ctime>
15 #include <chrono>
16 #include <algorithm>
17 
18 // mitk dependencies
19 #include "mitkUSDiPhASDevice.h"
21 #include <mitkIOUtil.h>
23 #include "ITKUltrasound/itkBModeImageFilter.h"
24 #include "mitkImageCast.h"
25 #include "mitkITKImageImport.h"
26 
27 // itk dependencies
28 #include "itkImage.h"
29 #include "itkResampleImageFilter.h"
30 #include "itkCastImageFilter.h"
31 #include "itkCropImageFilter.h"
32 #include "itkRescaleIntensityImageFilter.h"
33 #include "itkIntensityWindowingImageFilter.h"
34 #include <itkIndex.h>
35 #include "itkMultiplyImageFilter.h"
36 
38  : m_Device(device),
39  m_StartTime(((float)std::clock()) / CLOCKS_PER_SEC),
40  m_UseGUIOutPut(false),
41  m_DataType(DataType::Image_uChar),
42  m_GUIOutput(nullptr),
43  m_UseBModeFilter(false),
44  m_CurrentlyRecording(false),
45  m_DataTypeModified(true),
46  m_DataTypeNext(DataType::Image_uChar),
47  m_CurrentImageTimestamp(0),
48  m_PyroConnected(false),
49  m_ImageTimestampBuffer(),
50  m_VerticalSpacing(0),
51  m_UseBModeFilterModified(false),
52  m_UseBModeFilterNext(false),
53  m_ScatteringCoefficientModified(false),
54  m_CompensateForScatteringModified(false),
55  m_VerticalSpacingModified(false),
56  m_ScatteringCoefficient(15),
57  m_CompensateForScattering(false),
58  m_CompensateEnergy(false),
59  m_CompensateEnergyNext(false),
60  m_CompensateEnergyModified(false)
61 {
62  m_BufferSize = 100;
65  m_ImageBuffer.insert(m_ImageBuffer.begin(), m_BufferSize, nullptr);
66 
67  us::ModuleResource resourceFile;
68  std::string name;
70  for (int i = 5; i <= 25; ++i)
71  {
72  name = "c:\\HomogeneousScatteringAssumptions\\Scattering" + std::to_string(i) + ".nrrd";
73 
74  m_FluenceCompOriginal.push_back(mitk::IOUtil::Load<mitk::Image>(name));
75  }
76 
78  m_FluenceCompResizedItk.insert(m_FluenceCompResizedItk.begin(), 26, itk::Image<float,3>::New());
79 }
80 
82 {
83  // close the pyro
84  MITK_INFO("Pyro Debug") << "StopDataAcquisition: " << m_Pyro->StopDataAcquisition();
85  MITK_INFO("Pyro Debug") << "CloseConnection: " << m_Pyro->CloseConnection();
86  m_PyroConnected = false;
87  m_Pyro = nullptr;
88 }
89 
90 void mitk::USDiPhASImageSource::GetNextRawImage(std::vector<mitk::Image::Pointer>& imageVector)
91 {
92  // modify all settings that have been changed here, so we don't get multithreading issues
94  {
96  m_DataTypeModified = false;
98  }
100  {
102  m_UseBModeFilterModified = false;
103  }
105  {
108  }
110  {
113  }
115  {
118  }
120  {
123  }
124 
125  if (imageVector.size() != 2)
126  {
127  imageVector.resize(2);
128  }
129 
130  // make sure image is nullptr
131  mitk::Image::Pointer image = nullptr;
132  float ImageEnergyValue = 0;
133 
134  for (int i = 100; i > 90 && ImageEnergyValue <= 0; --i)
135  {
136  if (m_ImageTimestampBuffer[(m_LastWrittenImage + i) % 100] != 0)
137  {
138  ImageEnergyValue = m_Pyro->GetClosestEnergyInmJ(m_ImageTimestampBuffer[(m_LastWrittenImage + i) % 100]);
139  if (ImageEnergyValue > 0) {
140  image = &(*m_ImageBuffer[(m_LastWrittenImage + i) % 100]);
141  }
142  }
143  }
144  // if we did not get any usable Energy value, compensate using this default value
145  if (image == nullptr)
146  {
147  image = &(*m_ImageBuffer[m_LastWrittenImage]);
148  ImageEnergyValue = 40;
149  if (image == nullptr)
150  return;
151  }
152 
153  // do image processing before displaying it
154  if (image.IsNotNull())
155  {
156  itkFloatImageType::Pointer itkImage;
157 
158  mitk::CastToItkImage(image, itkImage);
159  image = mitk::GrabItkImageMemory(itkImage); //thereby using float images
160  image = CutOffTop(image, 165);
161 
162  // now apply filters to the image, if the options have been selected.
163  if ((m_CompensateForScattering || m_UseBModeFilter) && m_DataType == DataType::Beamformed_Short)
164  {
165  if (m_Device->GetScanMode().beamformingAlgorithm == Beamforming::PlaneWaveCompound) // this is for ultrasound only mode
166  {
167  if (m_UseBModeFilter)
168  {
169  image = ApplyBmodeFilter(image, true);
170  if (m_VerticalSpacing)
172  }
173  }
174 
175  else
176  {
177  Image::Pointer imagePA = Image::New();
178  unsigned int dim[] = { image->GetDimension(0),image->GetDimension(1),1};
179  imagePA->Initialize(image->GetPixelType(), 3, dim);
180  imagePA->SetGeometry(image->GetGeometry());
181 
182  Image::Pointer imageUS = Image::New();
183  imageUS->Initialize(image->GetPixelType(), 3, dim);
184  imageUS->SetGeometry(image->GetGeometry());
185 
186  ImageReadAccessor inputReadAccessorCopyPA(image, image->GetSliceData(0));
187  imagePA->SetSlice(inputReadAccessorCopyPA.GetData(), 0);
188  ImageReadAccessor inputReadAccessorCopyUS(image, image->GetSliceData(1));
189  imageUS->SetSlice(inputReadAccessorCopyUS.GetData(), 0);
190 
191  // first, seperate the PA image from the USImages
192 
193  // then, we compensate the PAImage using our ImageEnergyValue
195  imagePA = MultiplyImage(imagePA, 1/ImageEnergyValue); // TODO: add the correct prefactor here!!!!
196 
197  // now we apply the BModeFilter
198  if (m_UseBModeFilter)
199  {
200  imageUS = ApplyBmodeFilter(imageUS, true); // the US Images get a logarithmic filter
201  imagePA = ApplyBmodeFilter(imagePA, false);
202  }
203 
204  ImageReadAccessor inputReadAccessorPA(imagePA, imagePA->GetSliceData(0));
205  image->SetSlice(inputReadAccessorPA.GetData(), 0);
206  ImageReadAccessor inputReadAccessorUS(imageUS, imageUS->GetSliceData(0));
207  image->SetSlice(inputReadAccessorUS.GetData(), 1);
208  if (m_VerticalSpacing)
209  {
211  }
212 
213  // and finally the scattering corrections
215  {
216  auto curResizeImage = m_FluenceCompResized.at(m_ScatteringCoefficient); // just for convenience
217 
218  // update the fluence reference images!
219  bool doResampling = image->GetDimension(0) != curResizeImage->GetDimension(0) || image->GetDimension(1) != curResizeImage->GetDimension(1)
220  || image->GetGeometry()->GetSpacing()[0] != curResizeImage->GetGeometry()->GetSpacing()[0] || image->GetGeometry()->GetSpacing()[1] != curResizeImage->GetGeometry()->GetSpacing()[1];
221  if (doResampling)
222  {
223  curResizeImage = ApplyResampling(m_FluenceCompOriginal.at(m_ScatteringCoefficient), image->GetGeometry()->GetSpacing(), image->GetDimensions());
224 
225  double* rawOutputData = new double[image->GetDimension(0)*image->GetDimension(1)];
226  double* rawScatteringData = (double*)curResizeImage->GetData();
227  int sizeRawScatteringData = curResizeImage->GetDimension(0) * curResizeImage->GetDimension(1);
228  int imageSize = image->GetDimension(0)*image->GetDimension(1);
229 
230  //everything above 1.5mm is still inside the transducer; therefore the fluence compensation image has to be positioned a little lower
231  float upperCutoffmm = 1.5;
232  int lowerBound = std::round(upperCutoffmm / image->GetGeometry()->GetSpacing()[1])*image->GetDimension(0);
233  int upperBound = lowerBound + sizeRawScatteringData;
234 
235  for (int i = 0; i < lowerBound && i < imageSize; ++i)
236  {
237  rawOutputData[i] = 0; // everything than cannot be compensated shall be treated as garbage, here the upper 0.15mm
238  }
239  for (int i = lowerBound; i < upperBound && i < imageSize; ++i)
240  {
241  rawOutputData[i] = 1 / rawScatteringData[i-lowerBound];
242  }
243  for (int i = upperBound; i < imageSize; ++i)
244  {
245  rawOutputData[i] = 0; // everything than cannot be compensated shall be treated as garbage
246  }
247 
248 
249  unsigned int dim[] = { image->GetDimension(0), image->GetDimension(1), 1 };
250  curResizeImage->Initialize(mitk::MakeScalarPixelType<double>(), 3, dim);
251  curResizeImage->SetGeometry(image->GetGeometry());
252  curResizeImage->SetSlice(rawOutputData,0);
253 
254  delete[] rawOutputData;
255 
258 
259  MITK_INFO << "Resized a fluence image.";
260  }
261  // actually apply the scattering compensation
263  ImageReadAccessor inputReadAccessorPA(imagePA, imagePA->GetSliceData(0));
264  image->SetSlice(inputReadAccessorPA.GetData(), 0);
265  }
266  }
267  }
268 
269  //TODO: completely rewrite this mess
270 
271  imageVector[0] = Image::New();
272  unsigned int dim[] = { image->GetDimension(0),image->GetDimension(1),1 };
273  imageVector[0]->Initialize(image->GetPixelType(), 3, dim);
274  imageVector[0]->SetGeometry(image->GetGeometry());
275 
276  imageVector[1] = Image::New();
277  imageVector[1]->Initialize(image->GetPixelType(), 3, dim);
278  imageVector[1]->SetGeometry(image->GetGeometry());
279 
280  ImageReadAccessor inputReadAccessorCopyPA(image, image->GetSliceData(0));
281  imageVector[0]->SetSlice(inputReadAccessorCopyPA.GetData(), 0);
282  ImageReadAccessor inputReadAccessorCopyUS(image, image->GetSliceData(1));
283  imageVector[1]->SetSlice(inputReadAccessorCopyUS.GetData(), 0);
284  }
285 }
286 
288 {
289  // we use this seperate ApplyBmodeFilter Method for processing of two-dimensional images
290 
291  // the image needs to be of floating point type for the envelope filter to work; the casting is done automatically by the CastToItkImage
292 
294  BModeFilterType::Pointer bModeFilter = BModeFilterType::New(); // LogFilter
295 
297  PhotoacousticBModeImageFilter::Pointer photoacousticBModeFilter = PhotoacousticBModeImageFilter::New(); // No LogFilter
298 
299  itkFloatImageType::Pointer itkImage;
300  itkFloatImageType::Pointer bmode;
301  mitk::CastToItkImage(image, itkImage);
302 
303  if (useLogFilter)
304  {
305  bModeFilter->SetInput(itkImage);
306  bModeFilter->SetDirection(1);
307  bmode = bModeFilter->GetOutput();
308  }
309  else
310  {
311  photoacousticBModeFilter->SetInput(itkImage);
312  photoacousticBModeFilter->SetDirection(1);
313  bmode = photoacousticBModeFilter->GetOutput();
314  }
315  return mitk::GrabItkImageMemory(bmode);
316 }
317 
319 {
320  typedef itk::CropImageFilter < itkFloatImageType, itkFloatImageType > CutImageFilter;
321  itkFloatImageType::SizeType cropSize;
322  itkFloatImageType::Pointer itkImage;
323  mitk::CastToItkImage(image, itkImage);
324 
325  cropSize[0] = 0;
326  if(itkImage->GetLargestPossibleRegion().GetSize()[1] == 2048)
327  cropSize[1] = cutOffSize;
328  else
329  cropSize[1] = 0;
330  cropSize[2] = 0;
331  CutImageFilter::Pointer cutOffFilter = CutImageFilter::New();
332  cutOffFilter->SetInput(itkImage);
333  cutOffFilter->SetLowerBoundaryCropSize(cropSize);
334  cutOffFilter->UpdateLargestPossibleRegion();
335  return mitk::GrabItkImageMemory(cutOffFilter->GetOutput());
336 }
337 
339 {
340  typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter;
341  ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New();
342 
343  itkFloatImageType::Pointer itkImage;
344 
345  mitk::CastToItkImage(image, itkImage);
346  itkFloatImageType::SpacingType outputSpacing;
347  itkFloatImageType::SizeType inputSize = itkImage->GetLargestPossibleRegion().GetSize();
348  itkFloatImageType::SizeType outputSize = inputSize;
349 
350  outputSpacing[0] = itkImage->GetSpacing()[0] * (static_cast<double>(inputSize[0]) / static_cast<double>(outputSize[0]));
351  outputSpacing[1] = verticalSpacing;
352  outputSpacing[2] = itkImage->GetSpacing()[2];
353 
354  outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1];
355 
356  typedef itk::IdentityTransform<double, 3> TransformType;
357  resampleImageFilter->SetInput(itkImage);
358  resampleImageFilter->SetSize(outputSize);
359  resampleImageFilter->SetOutputSpacing(outputSpacing);
360  resampleImageFilter->SetTransform(TransformType::New());
361 
362  resampleImageFilter->UpdateLargestPossibleRegion();
363  return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput());
364 }
365 
367 {
368  typedef itk::MultiplyImageFilter <itkFloatImageType, itkFloatImageType > MultiplyImageFilterType;
369 
370  itkFloatImageType::Pointer itkImage;
371  mitk::CastToItkImage(inputImage, itkImage);
372 
373  MultiplyImageFilterType::Pointer multiplyFilter = MultiplyImageFilterType::New();
374  multiplyFilter->SetInput1(itkImage);
375  multiplyFilter->SetInput2(m_FluenceCompResizedItk.at(m_ScatteringCoefficient));
376 
377  return mitk::GrabItkImageMemory(multiplyFilter->GetOutput());
378 }
379 
381 {
382  typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter;
383  ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New();
384 
385  itkFloatImageType::Pointer itkImage;
386 
387  mitk::CastToItkImage(inputImage, itkImage);
388 
389  itkFloatImageType::SpacingType outputSpacingItk;
390  itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize();
391  itkFloatImageType::SizeType outputSizeItk = inputSizeItk;
392  itkFloatImageType::SpacingType inputSpacing = itkImage->GetSpacing();
393 
394  outputSizeItk[0] = outputSize[0];
395  outputSizeItk[1] = 10*(inputSpacing[1] * inputSizeItk[1]) / (outputSpacing[1]);
396  outputSizeItk[2] = 1;
397 
398  outputSpacingItk[0] = 0.996 * inputSpacing[0] * (static_cast<double>(inputSizeItk[0]) / static_cast<double>(outputSizeItk[0])); // TODO: find out why the spacing is not correct, so we need that factor; ?!?!
399  outputSpacingItk[1] = inputSpacing[1] * (static_cast<double>(inputSizeItk[1]) / static_cast<double>(outputSizeItk[1]));
400  outputSpacingItk[2] = outputSpacing[2];
401 
402  typedef itk::IdentityTransform<double, 3> TransformType;
403  resampleImageFilter->SetInput(itkImage);
404  resampleImageFilter->SetSize(outputSizeItk);
405  resampleImageFilter->SetOutputSpacing(outputSpacingItk);
406  resampleImageFilter->SetTransform(TransformType::New());
407 
408  resampleImageFilter->UpdateLargestPossibleRegion();
409  return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput());
410 }
411 
413 {
414  typedef itk::MultiplyImageFilter <itkFloatImageType, itkFloatImageType > MultiplyImageFilterType;
415 
416  itkFloatImageType::Pointer itkImage;
417  mitk::CastToItkImage(inputImage, itkImage);
418 
419  MultiplyImageFilterType::Pointer multiplyFilter = MultiplyImageFilterType::New();
420  multiplyFilter->SetInput1(itkImage);
421  multiplyFilter->SetConstant(value);
422 
423  return mitk::GrabItkImageMemory(multiplyFilter->GetOutput());
424 }
425 
427  short* rfDataChannelData,
428  int& channelDataChannelsPerDataset,
429  int& channelDataSamplesPerChannel,
430  int& channelDataTotalDatasets,
431 
432  short* rfDataArrayBeamformed,
433  int& beamformedLines,
434  int& beamformedSamples,
435  int& beamformedTotalDatasets,
436 
437  unsigned char* imageData,
438  int& imageWidth,
439  int& imageHeight,
440  int& imageBytesPerPixel,
441  int& imageSetsTotal,
442 
443  double& timeStamp)
444 {
445  if (m_DataTypeModified)
446  return;
447 
448  if (!m_PyroConnected)
449  {
451  MITK_INFO << "[Pyro Debug] OpenConnection: " << m_Pyro->OpenConnection();
452  MITK_INFO << "[Pyro Debug] StartDataAcquisition: " << m_Pyro->StartDataAcquisition();
453  m_PyroConnected = true;
454  }
455 
456  bool writeImage = ((m_DataType == DataType::Image_uChar) && (imageData != nullptr)) || ((m_DataType == DataType::Beamformed_Short) && (rfDataArrayBeamformed != nullptr));
457  if (writeImage)
458  {
459  //get the timestamp we might save later on
460  m_CurrentImageTimestamp = std::chrono::high_resolution_clock::now().time_since_epoch().count();
461 
462  // create a new image and initialize it
464 
465  switch (m_DataType)
466  {
467  case DataType::Image_uChar: {
468  m_ImageDimensions[0] = imageWidth;
469  m_ImageDimensions[1] = imageHeight;
470  m_ImageDimensions[2] = imageSetsTotal;
471  image->Initialize(mitk::MakeScalarPixelType<unsigned char>(), 3, m_ImageDimensions);
472  break;
473  }
474  case DataType::Beamformed_Short: {
475  m_ImageDimensions[0] = beamformedLines;
476  m_ImageDimensions[1] = beamformedSamples;
477  m_ImageDimensions[2] = beamformedTotalDatasets;
478  image->Initialize(mitk::MakeScalarPixelType<short>(), 3, m_ImageDimensions);
479  break;
480  }
481  }
482  image->GetGeometry()->SetSpacing(m_ImageSpacing);
483  image->GetGeometry()->Modified();
484 
485  // write the given buffer into the image
486  switch (m_DataType)
487  {
488  case DataType::Image_uChar: {
489  for (unsigned char i = 0; i < imageSetsTotal; i++) {
490  image->SetSlice(&imageData[i*imageHeight*imageWidth], i);
491  }
492  break;
493  }
494 
495  case DataType::Beamformed_Short: {
496  short* flipme = new short[beamformedLines*beamformedSamples*beamformedTotalDatasets];
497  int pixelsPerImage = beamformedLines*beamformedSamples;
498 
499  for (unsigned char currentSet = 0; currentSet < beamformedTotalDatasets; currentSet++)
500  {
501  for (unsigned short sample = 0; sample < beamformedSamples; sample++)
502  {
503  for (unsigned short line = 0; line < beamformedLines; line++)
504  {
505  flipme[sample*beamformedLines + line + pixelsPerImage*currentSet]
506  = rfDataArrayBeamformed[line*beamformedSamples + sample + pixelsPerImage*currentSet];
507  }
508  } // the beamformed pa image is flipped by 90 degrees; we need to flip it manually
509  }
510 
511  for (unsigned char i = 0; i < beamformedTotalDatasets; i++) {
512  image->SetSlice(&flipme[i*beamformedLines*beamformedSamples], i);
513  // set every image to a different slice
514  }
515 
516  delete[] flipme;
517  break;
518  }
519  }
520 
521  if (m_SavingSettings.saveRaw && m_CurrentlyRecording && rfDataChannelData != nullptr)
522  {
523  unsigned int dim[3];
524  dim[0] = channelDataChannelsPerDataset;
525  dim[1] = channelDataSamplesPerChannel;
526  dim[2] = 1;
527  short offset = m_Device->GetScanMode().accumulation * 2048;
528 
529  short* noOffset = new short[channelDataChannelsPerDataset*channelDataSamplesPerChannel*channelDataTotalDatasets];
530  for (unsigned char set = 0; set < 1; ++set)// channelDataTotalDatasets; ++set) // we ignore the raw US images for now
531  {
532  for (unsigned short sam = 0; sam < channelDataSamplesPerChannel; ++sam)
533  {
534  for (unsigned short chan = 0; chan < channelDataChannelsPerDataset; ++chan)
535  {
536  noOffset[set*channelDataSamplesPerChannel*channelDataChannelsPerDataset + sam * channelDataChannelsPerDataset + chan] =
537  rfDataChannelData[set*channelDataSamplesPerChannel*channelDataChannelsPerDataset + sam * channelDataChannelsPerDataset + chan] - offset; // this offset in the raw Images is given by the API...
538  }
539  }
540  }
541 
542  // save the raw images when recording
543  for (unsigned char i = 0; i < 1; ++i)// channelDataTotalDatasets; ++i) // we ignore the raw US images for now
544  {
546  rawImage->Initialize(mitk::MakeScalarPixelType<short>(), 3, dim);
547 
548  rawImage->SetSlice(&noOffset[i*channelDataChannelsPerDataset*channelDataSamplesPerChannel]);
549 
550  float& recordTime = m_Device->GetScanMode().receivePhaseLengthSeconds;
551  int& speedOfSound = m_Device->GetScanMode().averageSpeedOfSound;
552 
553  mitk::Vector3D rawSpacing;
554  rawSpacing[0] = m_Device->GetScanMode().transducerPitchMeter * 1000; // save in mm
555  rawSpacing[1] = recordTime / channelDataSamplesPerChannel * 1000000; // save in us
556  rawSpacing[2] = 1;
557 
558  rawImage->GetGeometry()->SetSpacing(rawSpacing);
559  rawImage->GetGeometry()->Modified();
560 
561  m_RawRecordedImages.push_back(rawImage);
562  }
563 
564  delete[] noOffset;
565  }
566 
567  itk::Index<3> pixel = { {
568  (itk::Index<3>::IndexValueType)(image->GetDimension(0) / 2),
569  (itk::Index<3>::IndexValueType)(22.0/532.0*m_Device->GetScanMode().reconstructionSamplesPerLine),
570  0 } }; //22/532*2048 = 84
571  if (!m_Pyro->IsSyncDelaySet() &&(image->GetPixelValueByIndex(pixel) < -30)) // #MagicNumber
572  {
573  MITK_INFO << "Setting SyncDelay now";
574  m_Pyro->SetSyncDelay(m_CurrentImageTimestamp);
575  }
576 
579  m_LastWrittenImage = (m_LastWrittenImage + 1) % m_BufferSize;
580 
581  // if the user decides to start recording, we feed the vector the generated images
582  if (m_CurrentlyRecording) {
583  for (unsigned char index = 0; index < image->GetDimension(2); ++index)
584  {
585  if (image->IsSliceSet(index))
586  {
587  m_RecordedImages.push_back(Image::New());
588  unsigned int dim[] = { image ->GetDimension(0), image->GetDimension(1), 1};
589  m_RecordedImages.back()->Initialize(image->GetPixelType(), 3, dim);
590  m_RecordedImages.back()->SetGeometry(image->GetGeometry());
591 
592  mitk::ImageReadAccessor inputReadAccessor(image, image->GetSliceData(index));
593  m_RecordedImages.back()->SetSlice(inputReadAccessor.GetData(),0);
594  }
595  }
597  // save timestamps for each laser image!
598  }
599  }
600 }
601 
603 {
604  MITK_INFO << "Retreaving Image Geometry Information for Spacing...";
605  float& recordTime = m_Device->GetScanMode().receivePhaseLengthSeconds;
606  int& speedOfSound = m_Device->GetScanMode().averageSpeedOfSound;
607  float& pitch = m_Device->GetScanMode().reconstructedLinePitchMmOrAngleDegree;
608  int& reconstructionLines = m_Device->GetScanMode().reconstructionLines;
609 
610  switch (m_DataType)
611  {
612  case DataType::Image_uChar : {
613  int& imageWidth = m_Device->GetScanMode().imageWidth;
614  int& imageHeight = m_Device->GetScanMode().imageHeight;
615  m_ImageSpacing[0] = pitch * reconstructionLines / imageWidth;
616  m_ImageSpacing[1] = recordTime * speedOfSound / 2 * 1000 / imageHeight;
617  break;
618  }
619  case DataType::Beamformed_Short : {
620  int& imageWidth = reconstructionLines;
621  int& imageHeight = m_Device->GetScanMode().reconstructionSamplesPerLine;
622  m_ImageSpacing[0] = pitch;
623  m_ImageSpacing[1] = recordTime * speedOfSound / 2 * 1000 / imageHeight;
624  break;
625  }
626  }
627  m_ImageSpacing[2] = 1;
628 
629  MITK_INFO << "Retreaving Image Geometry Information for Spacing " << m_ImageSpacing[0] << " ... " << m_ImageSpacing[1] << " ... " << m_ImageSpacing[2] << " ...[DONE]";
630 }
631 
633 {
634  m_DataTypeModified = true;
635  m_DataTypeNext = dataT;
636 }
637 
639 {
641  m_UseBModeFilterNext = isSet;
642 }
643 
645 {
648 }
649 
651 {
654 }
655 
657 {
658  m_CompensateEnergyNext = compensate;
660 }
661 
663 {
664  if (dataT != m_DataType)
665  {
666  m_DataType = dataT;
667  MITK_INFO << "Setting new DataType..." << dataT;
668  switch (m_DataType)
669  {
670  case DataType::Image_uChar :
671  MITK_INFO << "height: " << m_Device->GetScanMode().imageHeight << " width: " << m_Device->GetScanMode().imageWidth;
672  break;
673  case DataType::Beamformed_Short :
674  MITK_INFO << "samples: " << m_Device->GetScanMode().reconstructionSamplesPerLine << " lines: " << m_Device->GetScanMode().reconstructionLines;
675  break;
676  }
677  }
678 }
679 
680 void mitk::USDiPhASImageSource::SetGUIOutput(std::function<void(QString)> out)
681 {
683  m_StartTime = ((float)std::clock()) / CLOCKS_PER_SEC; //wait till the callback is available again
684  m_UseGUIOutPut = false;
685 }
686 
688 {
689  m_UseBModeFilter = isSet;
690 }
691 
693 {
696 }
697 
699 {
700  m_SavingSettings = settings;
701 }
702 
703 // this is just a little function to set the filenames below right
704 inline void replaceAll(std::string& str, const std::string& from, const std::string& to) {
705  if (from.empty())
706  return;
707  size_t start_pos = 0;
708  while ((start_pos = str.find(from, start_pos)) != std::string::npos) {
709  str.replace(start_pos, from.length(), to);
710  start_pos += to.length(); // In case 'to' contains 'from', like replacing 'x' with 'yx'
711  }
712 }
713 
715 {
716  // start the recording process
717  if (record)
718  {
719  m_RecordedImages.clear();
720  m_RawRecordedImages.clear(); // we make sure there are no leftovers
721  m_ImageTimestampRecord.clear(); // also for the timestamps
722  m_PixelValues.clear(); // aaaand for the pixel values
723 
725  {
726  m_Device->GetScanMode().transferChannelData = true;
728  // set the raw Data to be transfered
729  }
730 
731  // tell the callback to start recording images
732  m_CurrentlyRecording = true;
733  }
734  // save images, end recording, and clean up
735  else
736  {
737  m_CurrentlyRecording = false;
738 
739  m_Device->GetScanMode().transferChannelData = false; // make sure raw Channel Data is not transferred anymore!
741 
742  // get the time and date, put them into a nice string and create a folder for the images
743  time_t time = std::time(nullptr);
744  time_t* timeptr = &time;
745  std::string currentDate = std::ctime(timeptr);
746  replaceAll(currentDate, ":", "-");
747  currentDate.pop_back();
748  //std::string MakeFolder = "mkdir \"c:/DiPhASImageData/" + currentDate + "\"";
749  //system(MakeFolder.c_str());
750 
751  // initialize file paths and the images
752  Image::Pointer PAImage = Image::New();
754  std::string pathPA = "c:\\ImageData\\" + currentDate + "-" + "PAbeamformed" + ".nrrd";
755  std::string pathUS = "c:\\ImageData\\" + currentDate + "-" + "USImages" + ".nrrd";
756  std::string pathTS = "c:\\ImageData\\" + currentDate + "-" + "ts" + ".csv";
757  std::string pathS = "c:\\ImageData\\" + currentDate + "-" + "Settings" + ".txt";
758 
759  // idon't forget the raw Images (if chosen to be saved)
760  Image::Pointer PAImageRaw = Image::New();
761  Image::Pointer USImageRaw = Image::New();
762  std::string pathPARaw = "c:\\ImageData\\" + currentDate + "-" + "PAraw" + ".nrrd";
763  std::string pathUSRaw = "c:\\ImageData\\" + currentDate + "-" + "USImagesRaw" + ".nrrd";
764 
765  if (m_Device->GetScanMode().beamformingAlgorithm == (int)Beamforming::Interleaved_OA_US) // save a PAImage if we used interleaved mode
766  {
767  // first, save the data, so the pyro does not aquire more unneccessary timestamps
768  m_Pyro->SaveData();
769 
770  // now order the images and save them
771  // the beamformed ones...
773  {
774  OrderImagesInterleaved(PAImage, USImage, m_RecordedImages, false);
775  mitk::IOUtil::Save(USImage, pathUS);
776  mitk::IOUtil::Save(PAImage, pathPA);
777  }
778 
779  // ...and the raw images
781  {
782  OrderImagesInterleaved(PAImageRaw, USImageRaw, m_RawRecordedImages, true);
783  // mitk::IOUtil::Save(USImageRaw, pathUSRaw);
784  mitk::IOUtil::Save(PAImageRaw, pathPARaw);
785  }
786 
787  // read the pixelvalues of the enveloped images at this position
788 
789  itk::Index<3> pixel = { {
790  (itk::Index<3>::IndexValueType)(m_RecordedImages.at(0)->GetDimension(0) / 2),
791  (itk::Index<3>::IndexValueType)(22.0 / 532.0*m_Device->GetScanMode().reconstructionSamplesPerLine),
792  0 } }; //22/532*2048 = 84
793 
794  GetPixelValues(pixel, m_PixelValues); // write the Pixelvalues to m_PixelValues
795 
796  // save the timestamps!
797  ofstream timestampFile;
798 
799  timestampFile.open(pathTS);
800  timestampFile << ",timestamp,pixelvalue"; // write the header
801 
802  for (int index = 0; index < m_ImageTimestampRecord.size(); ++index)
803  {
804  timestampFile << "\n" << index << "," << m_ImageTimestampRecord.at(index) << "," << m_PixelValues.at(index);
805  }
806  timestampFile.close();
807 
808  //save the settings!
809 
810  ofstream settingsFile;
811 
812  settingsFile.open(pathS);
813  auto& sM = m_Device->GetScanMode();
814 
815  settingsFile << "[General Parameters]\n";
816  settingsFile << "Scan Depth [mm] = " << sM.receivePhaseLengthSeconds * sM.averageSpeedOfSound / 2 * 1000 << "\n";
817  settingsFile << "Speed of Sound [m/s] = " << sM.averageSpeedOfSound << "\n";
818  settingsFile << "Excitation Frequency [MHz] = " << sM.transducerFrequencyHz/1000000 << "\n";
819  settingsFile << "Voltage [V] = " << sM.voltageV << "\n";
820  settingsFile << "TGC min = " << (int)sM.tgcdB[0] << " max = " << (int)sM.tgcdB[7] << "\n";
821 
822  settingsFile << "[Beamforming Parameters]\n";
823  settingsFile << "Reconstructed Lines = " << sM.reconstructionLines << "\n";
824  settingsFile << "Samples per Line = " << sM.reconstructionSamplesPerLine << "\n";
825 
826  settingsFile.close();
827  }
828  else if (m_Device->GetScanMode().beamformingAlgorithm == (int)Beamforming::PlaneWaveCompound) // save no PAImage if we used US only mode
829  {
831  mitk::IOUtil::Save(USImage, pathUS);
832  }
833 
834  m_PixelValues.clear();
835  m_RawRecordedImages.clear(); // clean up the pixel values
836  m_RecordedImages.clear(); // clean up the images
837  m_ImageTimestampRecord.clear(); // clean up the timestamps
838  }
839 }
840 
841 void mitk::USDiPhASImageSource::GetPixelValues(itk::Index<3> pixel, std::vector<float>& values)
842 {
843  unsigned int events = 2;
844  for (int index = 0; index < m_RecordedImages.size(); index += events) // omit sound images
845  {
847  image = ApplyBmodeFilter(image);
848  values.push_back(image.GetPointer()->GetPixelValueByIndex(pixel));
849  }
850 }
851 
852 void mitk::USDiPhASImageSource::OrderImagesInterleaved(Image::Pointer PAImage, Image::Pointer USImage, std::vector<Image::Pointer> recordedList, bool raw)
853 {
854  unsigned int width = 32;
855  unsigned int height = 32;
856  unsigned int events = m_Device->GetScanMode().transmitEventsCount + 1; // the PA event is not included in the transmitEvents, so we add 1 here
857  if (!raw)
858  events = 2; // the beamformed image array contains only the resulting image of multiple events
859 
860  if (raw)
861  {
862  width = recordedList.at(0)->GetDimension(0);
863  height = recordedList.at(0)->GetDimension(1);
864  }
865  else if (m_DataType == DataType::Beamformed_Short)
866  {
867  width = m_Device->GetScanMode().reconstructionLines;
868  height = m_Device->GetScanMode().reconstructionSamplesPerLine;
869  }
870  else if (m_DataType == DataType::Image_uChar)
871  {
872  width = m_Device->GetScanMode().imageWidth;
873  height = m_Device->GetScanMode().imageHeight;
874  }
875 
876  unsigned int dimLaser[] = { (unsigned int)width, (unsigned int)height, (unsigned int)(recordedList.size() / events)};
877  unsigned int dimSound[] = { (unsigned int)width, (unsigned int)height, (unsigned int)(recordedList.size() / events * (events-1))};
878 
879  PAImage->Initialize(recordedList.back()->GetPixelType(), 3, dimLaser);
880  PAImage->SetGeometry(recordedList.back()->GetGeometry());
881  USImage->Initialize(recordedList.back()->GetPixelType(), 3, dimSound);
882  USImage->SetGeometry(recordedList.back()->GetGeometry());
883 
884  for (int index = 0; index < recordedList.size(); ++index)
885  {
886  mitk::ImageReadAccessor inputReadAccessor(recordedList.at(index));
887 
888  if (index % events == 0)
889  {
890  PAImage->SetSlice(inputReadAccessor.GetData(), index / events);
891  }
892  else
893  {
894  if(!raw)
895  USImage->SetSlice(inputReadAccessor.GetData(), ((index - (index % events)) / events) + (index % events)-1);
896  }
897  }
898 }
899 
900 void mitk::USDiPhASImageSource::OrderImagesUltrasound(Image::Pointer USImage, std::vector<Image::Pointer> recordedList)
901 {
902  unsigned int width = 32;
903  unsigned int height = 32;
904  unsigned int events = m_Device->GetScanMode().transmitEventsCount;
905 
906  if (m_DataType == DataType::Beamformed_Short)
907  {
908  width = (unsigned int)m_Device->GetScanMode().reconstructionLines;
909  height = (unsigned int)m_Device->GetScanMode().reconstructionSamplesPerLine;
910  }
911  else if (m_DataType == DataType::Image_uChar)
912  {
913  width = (unsigned int)m_Device->GetScanMode().imageWidth;
914  height = (unsigned int)m_Device->GetScanMode().imageHeight;
915  }
916 
917  unsigned int dimSound[] = { (unsigned int)width, (unsigned int)height, (unsigned int)recordedList.size()};
918 
919  USImage->Initialize(recordedList.back()->GetPixelType(), 3, dimSound);
920  USImage->SetGeometry(recordedList.back()->GetGeometry());
921 
922  for (int index = 0; index < recordedList.size(); ++index)
923  {
924  mitk::ImageReadAccessor inputReadAccessor(recordedList.at(index));
925  USImage->SetSlice(inputReadAccessor.GetData(), index);
926  }
927 }
void ImageDataCallback(short *rfDataChannelData, int &channelDataChannelsPerDataset, int &channelDataSamplesPerChannel, int &channelDataTotalDatasets, short *rfDataArrayBeamformed, int &beamformedLines, int &beamformedSamples, int &beamformedTotalDatasets, unsigned char *imageData, int &imageWidth, int &imageHeight, int &imagePixelFormat, int &imageSetsTotal, double &timeStamp)
static char * line
Definition: svm.cpp:2870
std::vector< mitk::Image::Pointer > m_ImageBuffer
#define MITK_INFO
Definition: mitkLogMacros.h:18
void replaceAll(std::string &str, const std::string &from, const std::string &to)
void OrderImagesUltrasound(Image::Pointer USImage, std::vector< Image::Pointer > recordedList)
static Pointer New()
void SetGUIOutput(std::function< void(QString)> out)
mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, mitk::Vector3D outputSpacing, unsigned int outputSize[3])
mitk::USDiPhASDevice * m_Device
STL namespace.
std::vector< Image::Pointer > m_FluenceCompOriginal
void ModifyEnergyCompensation(bool compensate)
mitk::OphirPyro::Pointer m_Pyro
Image::Pointer GrabItkImageMemory(itk::SmartPointer< ItkOutputImageType > &itkimage, mitk::Image *mitkImage=nullptr, const BaseGeometry *geometry=nullptr, bool update=true)
Grabs the memory of an itk::Image (with a specific type) and puts it into an mitk::Image.The memory is managed by the mitk::Image after calling this function. The itk::Image remains valid until the mitk::Image decides to free the memory.
std::vector< float > m_PixelValues
This specialization of mitk::Image only appends necessary Metadata to an MITK image. Otherwise it can safely be treated like it&#39;s mother class. To generate an USImage from a standard mitkImage, call the appropriate constructor USImage(image::Pointer)
Definition: mitkUSImage.h:31
static Vector3D offset
std::vector< mitk::Image::Pointer > m_RawRecordedImages
std::vector< long long > m_ImageTimestampRecord
mitk::Image::Pointer ResampleOutputVertical(mitk::Image::Pointer image, float verticalSpacing=0.1)
mitk::Image::Pointer MultiplyImage(mitk::Image::Pointer inputImage, double value)
std::vector< long long > m_ImageTimestampBuffer
std::vector< Image::Pointer > m_FluenceCompResized
void ModifyCompensateForScattering(bool useIt)
void SetSavingSettings(SavingSettings settings)
std::function< void(QString)> m_GUIOutput
mitk::Image::Pointer image
static Pointer New()
mitk::Image::Pointer ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scatteringCoefficient)
std::vector< itk::Image< float, 3 >::Pointer > m_FluenceCompResizedItk
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
mitk::Image::Pointer CutOffTop(mitk::Image::Pointer image, int cutOffSize=165)
static void Save(const mitk::BaseData *data, const std::string &path, bool setPathProperty=false)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:774
void GetPixelValues(itk::Index< 3 > pixel, std::vector< float > &values)
ScanModeNative & GetScanMode()
virtual void GetNextRawImage(std::vector< mitk::Image::Pointer > &) override
Create an Photoacoustic B-Mode (Brightness-Mode) image from raw "RF" data. The RF&#39;s envelope is calcu...
ImageReadAccessor class to get locked read access for a particular image part.
USDiPhASImageSource(mitk::USDiPhASDevice *device)
std::vector< mitk::Image::Pointer > m_RecordedImages
Create an ultrasound B-Mode (Brightness-Mode) image from raw "RF" data. The RF&#39;s envelope is calculat...
mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer image, bool useLogFilter=false)
void OrderImagesInterleaved(Image::Pointer PAImage, Image::Pointer USImage, std::vector< Image::Pointer > recordedList, bool raw)