diff --git a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp
index 81f0520417..f030cc339e 100644
--- a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp
+++ b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp
@@ -1,232 +1,233 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkCommon.h>
 #include <chrono>
 #include <mitkIOUtil.h>
 #include <mitkCommandLineParser.h>
 #include <mitkException.h>
 
 #include <mitkPhotoacousticFilterService.h>
 #include <mitkBeamformingSettings.h>
 #include <mitkCastToFloatImageFilter.h>
 
 #include <itksys/SystemTools.hxx>
 
 struct InputParameters
 {
   mitk::Image::Pointer inputImage;
   std::string outputFilename;
   bool verbose;
   float speedOfSound;
   unsigned int cutoff;
   float angle;
   unsigned int samples;
   mitk::BeamformingSettings::BeamformingAlgorithm algorithm;
 };
 
 InputParameters parseInput(int argc, char* argv[])
 {
   mitkCommandLineParser parser;
   parser.setCategory("MITK-Photoacoustics");
   parser.setTitle("Mitk Photoacoustics Beamforming Tool");
   parser.setDescription("Reads a nrrd file as an input and applies a beamforming method as set with the parameters.");
   parser.setContributor("Computer Assisted Medical Interventions, DKFZ");
 
   parser.setArgumentPrefix("--", "-");
 
   parser.beginGroup("Required parameters");
   parser.addArgument(
     "inputImage", "i", mitkCommandLineParser::InputImage,
     "Input image (mitk::Image)", "input image (.nrrd file)",
     us::Any(), false);
   parser.addArgument(
     "output", "o", mitkCommandLineParser::OutputFile,
     "Output filename", "output image (.nrrd file)",
     us::Any(), false);
   parser.endGroup();
 
   parser.beginGroup("Optional parameters");
   parser.addArgument(
     "verbose", "v", mitkCommandLineParser::Bool,
     "Verbose Output", "Whether to produce verbose, or rather debug output. (default: false)");
   parser.addArgument(
     "speed-of-sound", "sos", mitkCommandLineParser::Float,
     "Speed of Sound [m/s]", "The average speed of sound as assumed for the reconstruction in [m/s]. (default: 1500)");
   parser.addArgument(
     "cutoff", "co", mitkCommandLineParser::Int,
     "cutoff margin on the top of the image [pixels]", "The number of pixels to be ignored for this filter in [pixels] (default: 0).");
   parser.addArgument(
     "angle", "a", mitkCommandLineParser::Float,
     "field of view of the transducer elements [degrees]", "The field of view of each individual transducer element [degrees] (default: 27).");
   parser.addArgument(
     "samples", "s", mitkCommandLineParser::Int,
     "samples per reconstruction line [pixels]", "The pixels along the y axis in the beamformed image [pixels] (default: 2048).");
   parser.addArgument(
     "algorithm", "alg", mitkCommandLineParser::String,
     "one of [\"DAS\", \"DMAS\", \"sDMAS\"]", "The beamforming algorithm to be used for reconstruction (default: DAS).");
   parser.endGroup();
 
   InputParameters input;
 
   std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
   if (parsedArgs.size() == 0)
     exit(-1);
 
   if (parsedArgs.count("verbose"))
   {
     input.verbose = true;
   }
   else
   {
     input.verbose = false;
   }
   MITK_INFO(input.verbose) << "### VERBOSE OUTPUT ENABLED ###";
 
   if (parsedArgs.count("inputImage"))
   {
     MITK_INFO(input.verbose) << "Reading input image...";
     input.inputImage = mitk::IOUtil::Load<mitk::Image>(us::any_cast<std::string>(parsedArgs["inputImage"]));
     MITK_INFO(input.verbose) << "Reading input image...[Done]";
   }
   else
   {
     mitkThrow() << "No input image given.";
   }
 
   if (parsedArgs.count("output"))
   {
     input.outputFilename = us::any_cast<std::string>(parsedArgs["output"]);
   }
   else
   {
     mitkThrow() << "No output image path given..";
   }
 
   if (parsedArgs.count("speed-of-sound"))
   {
     input.speedOfSound = us::any_cast<int>(parsedArgs["speed-of-sound"]);
   }
   else
   {
     input.speedOfSound = 1500;
   }
 
   if (parsedArgs.count("cutoff"))
   {
     input.cutoff = us::any_cast<int>(parsedArgs["cutoff"]);
   }
   else
   {
     input.cutoff = 0;
   }
 
   if (parsedArgs.count("angle"))
   {
     input.angle = us::any_cast<float>(parsedArgs["angle"]);
   }
   else
   {
     input.angle = 27;
   }
 
   if (parsedArgs.count("samples"))
   {
     input.samples = us::any_cast<int>(parsedArgs["samples"]);
   }
   else
   {
     input.samples = 2048;
   }
 
   if (parsedArgs.count("algorithm"))
   {
     std::string algorithm = us::any_cast<std::string>(parsedArgs["algorithm"]);
     MITK_INFO(input.verbose) << "Parsing algorithm: " << algorithm;
     if (algorithm == "DAS")
       input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS;
     else if (algorithm == "DMAS")
       input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS;
     else if (algorithm == "sDMAS")
       input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS;
     else
     {
       MITK_INFO(input.verbose) << "Not a valid beamforming algorithm: " << algorithm << " Reverting to DAS";
       input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS;
     }
 
     MITK_INFO(input.verbose) << "Sucessfully set algorithm: " << algorithm;
   }
   else
   {
     input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS;
     MITK_INFO(input.verbose) << "No matching algorithm found. Using DAS.";
   }
 
   return input;
 }
 
 mitk::BeamformingSettings::Pointer ParseSettings(InputParameters &input)
 {
   mitk::BeamformingSettings::Pointer outputSettings = mitk::BeamformingSettings::New(
     (float)(input.inputImage->GetGeometry()->GetSpacing()[0] / 1000),
     (float)(input.speedOfSound),
     (float)(input.inputImage->GetGeometry()->GetSpacing()[1] / 1000000),
     input.angle,
     true,
+    input.inputImage->GetDimension(1),
     input.inputImage->GetDimension(0),
     input.cutoff,
     false,
     (unsigned int*) nullptr,
     input.inputImage->GetDimensions(),
     false,
     mitk::BeamformingSettings::DelayCalc::Spherical,
     mitk::BeamformingSettings::Apodization::Box,
     input.inputImage->GetDimension(0),
     input.algorithm,
     false,
     0.0f,
     0.0f);
 
   return outputSettings;
 }
 
 int main(int argc, char * argv[])
 {
   auto input = parseInput(argc, argv);
 
   MITK_INFO(input.verbose) << "Beamforming input image...";
 
   mitk::PhotoacousticFilterService::Pointer m_BeamformingService = mitk::PhotoacousticFilterService::New();
 
   mitk::BeamformingSettings::Pointer settings = ParseSettings(input);
 
   mitk::CastToFloatImageFilter::Pointer castFilter = mitk::CastToFloatImageFilter::New();
   castFilter->SetInput(input.inputImage);
   castFilter->Update();
   auto floatImage = castFilter->GetOutput();
 
   auto output = m_BeamformingService->ApplyBeamforming(floatImage, settings);
   MITK_INFO(input.verbose) << "Applying BModeFilter to image...";
-  //output = m_BeamformingService->ApplyBmodeFilter(output, mitk::PhotoacousticFilterService::Abs, false, false, 0.3);
+  auto output2 = m_BeamformingService->ApplyBmodeFilter(output, mitk::PhotoacousticFilterService::Abs, false, false, 0.3);
   MITK_INFO(input.verbose) << "Applying BModeFilter to image...[Done]";
 
   MITK_INFO(input.verbose) << "Saving image...";
-  mitk::IOUtil::Save(output, input.outputFilename);
+  mitk::IOUtil::Save(output2, input.outputFilename);
   MITK_INFO(input.verbose) << "Saving image...[Done]";
 
   MITK_INFO(input.verbose) << "Beamforming input image...[Done]";
 }
diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp
index f9db6a7105..65164364df 100644
--- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp
+++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp
@@ -1,221 +1,222 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h"
 #include "usServiceReference.h"
 #include <mitkImageReadAccessor.h>
 
 #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN
 
 mitk::PhotoacousticOCLBModeFilter::PhotoacousticOCLBModeFilter()
   : m_PixelCalculation(NULL)
 {
   this->AddSourceFile("BModeAbs.cl");
   this->AddSourceFile("BModeAbsLog.cl");
 
   this->m_FilterID = "BModeFilter";
 
   this->Initialize();
 }
 
 mitk::PhotoacousticOCLBModeFilter::~PhotoacousticOCLBModeFilter()
 {
   if (this->m_PixelCalculation)
   {
     clReleaseKernel(m_PixelCalculation);
   }
 }
 
 void mitk::PhotoacousticOCLBModeFilter::Update()
 {
   //Check if context & program available
   if (!this->Initialize())
   {
     us::ServiceReference<OclResourceService> ref = GetModuleContext()->GetServiceReference<OclResourceService>();
     OclResourceService* resources = GetModuleContext()->GetService<OclResourceService>(ref);
 
     // clean-up also the resources
     resources->InvalidateStorage();
     mitkThrow() << "Filter is not initialized. Cannot update.";
   }
   else {
     // Execute
     this->Execute();
   }
 }
 
 void mitk::PhotoacousticOCLBModeFilter::Execute()
 {
   try
   {
     size_t outputSize = m_InputDim[0] * m_InputDim[1] * m_InputDim[2];
     this->InitExec(this->m_PixelCalculation, m_InputDim, outputSize, sizeof(float));
   }
   catch (const mitk::Exception& e)
   {
     MITK_ERROR << "Catched exception while initializing filter: " << e.what();
     return;
   }
 
   cl_int clErr;
   clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Size));
 
   CHECK_OCL_ERR(clErr);
 
   // execute the filter on a 3D NDRange
   this->ExecuteKernel(m_PixelCalculation, 3);
 
   // signalize the GPU-side data changed
   m_Output->Modified(GPU_DATA);
 }
 
 us::Module *mitk::PhotoacousticOCLBModeFilter::GetModule()
 {
   return us::GetModuleContext()->GetModule();
 }
 
 bool mitk::PhotoacousticOCLBModeFilter::Initialize()
 {
   bool buildErr = true;
   cl_int clErr = 0;
 
   if (OclFilter::Initialize())
   {
     if(m_UseLogFilter)
       this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckBmodeAbsLog", &clErr);
     else
       this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckBmodeAbs", &clErr);
     buildErr |= CHECK_OCL_ERR(clErr);
   }
   return (OclFilter::IsInitialized() && buildErr);
 }
 
 void mitk::PhotoacousticOCLBModeFilter::SetInput(mitk::Image::Pointer image)
 {
   OclDataSetToDataSetFilter::SetInput(image);
 
   m_InputImage = image;
   m_InputDim[0] = m_InputImage->GetDimension(0);
   m_InputDim[1] = m_InputImage->GetDimension(1);
   m_InputDim[2] = m_InputImage->GetDimension(2);
   m_Size = m_InputDim[0] * m_InputDim[1] * m_InputDim[2];
 }
 
 mitk::Image::Pointer mitk::PhotoacousticOCLBModeFilter::GetOutput()
 {
   mitk::Image::Pointer outputImage = mitk::Image::New();
 
   if (m_Output->IsModified(GPU_DATA))
   {
     void* pData = m_Output->TransferDataToCPU(m_CommandQue);
 
     const unsigned int dimension = 3;
     unsigned int dimensions[3] = { m_InputDim[0], m_InputDim[1], m_InputDim[2] };
 
     const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry();
 
     MITK_DEBUG << "Creating new MITK Image.";
 
     outputImage->Initialize(this->GetOutputType(), dimension, dimensions);
     outputImage->SetSpacing(p_slg->GetSpacing());
     outputImage->SetGeometry(m_InputImage->GetGeometry());
-    outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ReferenceMemory);
+    outputImage->SetImportVolume(pData, 0, 0, mitk::Image::CopyMemory);
+    delete[] pData;
   }
 
   MITK_DEBUG << "Image Initialized.";
 
   return outputImage;
 }
 
 #endif
 
 
 mitk::PhotoacousticBModeFilter::PhotoacousticBModeFilter() : m_UseLogFilter(false)
 {
   this->SetNumberOfIndexedInputs(1);
   this->SetNumberOfRequiredInputs(1);
 }
 
 mitk::PhotoacousticBModeFilter::~PhotoacousticBModeFilter()
 {
 }
 
 void mitk::PhotoacousticBModeFilter::GenerateInputRequestedRegion()
 {
   Superclass::GenerateInputRequestedRegion();
 
   mitk::Image* output = this->GetOutput();
   mitk::Image* input = const_cast<mitk::Image *> (this->GetInput());
   if (!output->IsInitialized())
   {
     return;
   }
 
   input->SetRequestedRegionToLargestPossibleRegion();
 
   //GenerateTimeInInputRegion(output, input);
 }
 
 void mitk::PhotoacousticBModeFilter::GenerateOutputInformation()
 {
   mitk::Image::ConstPointer input = this->GetInput();
   mitk::Image::Pointer output = this->GetOutput();
 
   if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
     return;
 
   itkDebugMacro(<< "GenerateOutputInformation()");
 
-  output->Initialize(input);
+  output->Initialize(input->GetPixelType(), input->GetDimension(), input->GetDimensions());
   output->GetGeometry()->SetSpacing(input->GetGeometry()->GetSpacing());
   output->GetGeometry()->Modified();
   output->SetPropertyList(input->GetPropertyList()->Clone());
 
   m_TimeOfHeaderInitialization.Modified();
 }
 
 void mitk::PhotoacousticBModeFilter::GenerateData()
 {
   GenerateOutputInformation();
   mitk::Image::Pointer input = this->GetInput();
   mitk::Image::Pointer output = this->GetOutput();
 
   if (!output->IsInitialized())
     return;
 
   mitk::ImageReadAccessor reader(input);
 
   unsigned int size = output->GetDimension(0) * output->GetDimension(1) * output->GetDimension(2);
 
-  float* InputData = (float*)const_cast<void*>(reader.GetData());
+  const float* InputData = (const float*)(reader.GetData());
   float* OutputData = new float[size];
   if(!m_UseLogFilter)
     for (unsigned int i = 0; i < size; ++i)
     {
       OutputData[i] = abs(InputData[i]);
     }
   else
   {
     for (unsigned int i = 0; i < size; ++i)
     {
       OutputData[i] = log(abs(InputData[i]));
     }
   }
 
-  output->SetImportVolume(OutputData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory);
-
+  output->SetImportVolume(OutputData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory);
+  delete[] OutputData;
   m_TimeOfHeaderInitialization.Modified();
 }
\ No newline at end of file
diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp
index ffa3db3902..ed8914d7c4 100644
--- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp
+++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp
@@ -1,441 +1,443 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPhotoacousticFilterService.h"
 #include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h"
 #include "../ITKFilter/itkPhotoacousticBModeImageFilter.h"
 #include "mitkImageCast.h"
 #include "mitkITKImageImport.h"
 #include "mitkBeamformingFilter.h"
 #include <chrono>
 #include <mitkIOUtil.h>
 #include <mitkCropImageFilter.h>
 #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h"
 #include "mitkConvert2Dto3DImageFilter.h"
 #include <mitkCastToFloatImageFilter.h>
 
 // itk dependencies
 #include "itkImage.h"
 #include "itkResampleImageFilter.h"
 #include "itkCastImageFilter.h"
 #include "itkCropImageFilter.h"
 #include "itkRescaleIntensityImageFilter.h"
 #include "itkIntensityWindowingImageFilter.h"
 #include <itkIndex.h>
 #include "itkMultiplyImageFilter.h"
 #include "itkBSplineInterpolateImageFunction.h"
 #include <mitkImageToItk.h>
 
 // needed itk image filters
 #include "mitkITKImageImport.h"
 #include "itkFFTShiftImageFilter.h"
 #include "itkMultiplyImageFilter.h"
 #include "itkComplexToModulusImageFilter.h"
 #include <itkAddImageFilter.h>
 #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h"
 #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h"
 
 mitk::PhotoacousticFilterService::PhotoacousticFilterService()
 {
   MITK_INFO << "[PhotoacousticFilterService] created filter service";
 }
 
 mitk::PhotoacousticFilterService::~PhotoacousticFilterService()
 {
   MITK_INFO << "[PhotoacousticFilterService] destructed filter service";
 }
 
 mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseGPU, bool UseLogFilter, float resampleSpacing)
 {
   // the image needs to be of floating point type for the envelope filter to work; the casting is done automatically by the CastToItkImage
   typedef itk::Image< float, 3 > itkFloatImageType;
   typedef itk::IdentityTransform<double, 3> TransformType;
 
+  if (inputImage.IsNull() || !(inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)"))
+  {
+    MITK_ERROR << "BMode Filter can only handle float image types.";
+    mitkThrow() << "BMode Filter can only handle float image types.";
+  }
+
   if (method == BModeMethod::Abs)
   {
-    mitk::Image::Pointer input;
+    mitk::Image::Pointer input = inputImage;
     mitk::Image::Pointer out;
-    if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)")
-      input = inputImage;
-    else
-      input = ApplyCropping(inputImage, 0, 0, 0, 0, 0, 0);
 
     if (!UseGPU)
     {
       PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New();
       filter->SetParameters(UseLogFilter);
       filter->SetInput(input);
       filter->Update();
 
       out = filter->GetOutput();
 
       if (resampleSpacing == 0)
         return out;
     }
 #ifdef PHOTOACOUSTICS_USE_GPU
     else
     {
       PhotoacousticOCLBModeFilter::Pointer filter = PhotoacousticOCLBModeFilter::New();
       filter->SetParameters(UseLogFilter);
       filter->SetInput(input);
       filter->Update();
 
       out = filter->GetOutput();
 
       if (resampleSpacing == 0)
         return out;
     }
 #endif
 
     typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter;
     ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New();
 
-    itkFloatImageType::Pointer itkImage;
+    itkFloatImageType::Pointer itkImage = itkFloatImageType::New();
     mitk::CastToItkImage(out, itkImage);
     itkFloatImageType::SpacingType outputSpacing;
     itkFloatImageType::SizeType inputSize = itkImage->GetLargestPossibleRegion().GetSize();
     itkFloatImageType::SizeType outputSize = inputSize;
 
     outputSpacing[0] = itkImage->GetSpacing()[0];
     outputSpacing[1] = resampleSpacing;
     outputSpacing[2] = itkImage->GetSpacing()[2];
 
     outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1];
 
     typedef itk::IdentityTransform<double, 3> TransformType;
     resampleImageFilter->SetInput(itkImage);
     resampleImageFilter->SetSize(outputSize);
     resampleImageFilter->SetOutputSpacing(outputSpacing);
     resampleImageFilter->SetTransform(TransformType::New());
 
     resampleImageFilter->UpdateLargestPossibleRegion();
     return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput());
   }
   else if (method == BModeMethod::EnvelopeDetection)
   {
     typedef itk::BModeImageFilter < itkFloatImageType, itkFloatImageType > BModeFilterType;
     BModeFilterType::Pointer bModeFilter = BModeFilterType::New();  // LogFilter
 
     typedef itk::PhotoacousticBModeImageFilter < itkFloatImageType, itkFloatImageType > PhotoacousticBModeImageFilter;
     PhotoacousticBModeImageFilter::Pointer photoacousticBModeFilter = PhotoacousticBModeImageFilter::New(); // No LogFilter
 
     typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter;
     ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New();
 
     itkFloatImageType::Pointer itkImage;
 
     mitk::CastToItkImage(inputImage, itkImage);
 
     itkFloatImageType::Pointer bmode;
 
     if (UseLogFilter)
     {
       bModeFilter->SetInput(itkImage);
       bModeFilter->SetDirection(1);
       bmode = bModeFilter->GetOutput();
     }
     else
     {
       photoacousticBModeFilter->SetInput(itkImage);
       photoacousticBModeFilter->SetDirection(1);
       bmode = photoacousticBModeFilter->GetOutput();
     }
 
     // resampleSpacing == 0 means: do no resampling
     if (resampleSpacing == 0)
     {
       return mitk::GrabItkImageMemory(bmode);
     }
 
     itkFloatImageType::SpacingType outputSpacing;
     itkFloatImageType::SizeType inputSize = itkImage->GetLargestPossibleRegion().GetSize();
     itkFloatImageType::SizeType outputSize = inputSize;
 
     outputSpacing[0] = itkImage->GetSpacing()[0];
     outputSpacing[1] = resampleSpacing;
     outputSpacing[2] = itkImage->GetSpacing()[2];
 
     outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1];
 
     resampleImageFilter->SetInput(bmode);
     resampleImageFilter->SetSize(outputSize);
     resampleImageFilter->SetOutputSpacing(outputSpacing);
     resampleImageFilter->SetTransform(TransformType::New());
 
     resampleImageFilter->UpdateLargestPossibleRegion();
     return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput());
   }
   return nullptr;
 }
 
 mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2])
 {
   typedef itk::Image< float, 3 > itkFloatImageType;
 
   typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter;
   ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New();
   typedef itk::LinearInterpolateImageFunction<itkFloatImageType, double> T_Interpolator;
   itkFloatImageType::Pointer itkImage;
 
   mitk::CastToItkImage(inputImage, itkImage);
 
   itkFloatImageType::SpacingType outputSpacingItk;
   itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize();
   itkFloatImageType::SizeType outputSizeItk = inputSizeItk;
 
   outputSizeItk[0] = outputSize[0];
   outputSizeItk[1] = outputSize[1];
   outputSizeItk[2] = inputSizeItk[2];
 
   outputSpacingItk[0] = itkImage->GetSpacing()[0] * (static_cast<double>(inputSizeItk[0]) / static_cast<double>(outputSizeItk[0]));
   outputSpacingItk[1] = itkImage->GetSpacing()[1] * (static_cast<double>(inputSizeItk[1]) / static_cast<double>(outputSizeItk[1]));
   outputSpacingItk[2] = itkImage->GetSpacing()[2];
 
   typedef itk::IdentityTransform<double, 3> TransformType;
   T_Interpolator::Pointer _pInterpolator = T_Interpolator::New();
   resampleImageFilter->SetInput(itkImage);
   resampleImageFilter->SetSize(outputSizeItk);
   resampleImageFilter->SetOutputSpacing(outputSpacingItk);
   resampleImageFilter->SetTransform(TransformType::New());
   resampleImageFilter->SetInterpolator(_pInterpolator);
 
   resampleImageFilter->UpdateLargestPossibleRegion();
   return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput());
 }
 
 mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyCropping(mitk::Image::Pointer inputImage,
   int above, int below,
   int right, int left,
   int zStart, int zEnd)
 {
   mitk::CropImageFilter::Pointer cropImageFilter = mitk::CropImageFilter::New();
   cropImageFilter->SetInput(inputImage);
   cropImageFilter->SetXPixelsCropStart(left);
   cropImageFilter->SetXPixelsCropEnd(right);
   cropImageFilter->SetYPixelsCropStart(above);
   cropImageFilter->SetYPixelsCropEnd(below);
   cropImageFilter->SetZPixelsCropStart(zStart);
   cropImageFilter->SetZPixelsCropEnd(zEnd);
   cropImageFilter->Update();
   return cropImageFilter->GetOutput();
 }
 
 mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBeamforming(mitk::Image::Pointer inputImage,
   BeamformingSettings::Pointer config, std::function<void(int, std::string)> progressHandle)
 {
   Image::Pointer processedImage = mitk::Image::New();
 
   if (inputImage->GetDimension() != 3)
   {
     mitk::Convert2Dto3DImageFilter::Pointer dimensionImageFilter = mitk::Convert2Dto3DImageFilter::New();
     dimensionImageFilter->SetInput(inputImage);
     dimensionImageFilter->Update();
     processedImage = dimensionImageFilter->GetOutput();
   }
   else
   {
     processedImage = inputImage;
   }
 
   // perform the beamforming
   m_BeamformingFilter = mitk::BeamformingFilter::New(config);
   m_BeamformingFilter->SetInput(processedImage);
   m_BeamformingFilter->SetProgressHandle(progressHandle);
   m_BeamformingFilter->UpdateLargestPossibleRegion();
 
   processedImage = m_BeamformingFilter->GetOutput();
 
   return processedImage;
 }
 
 mitk::Image::Pointer mitk::PhotoacousticFilterService::BandpassFilter(mitk::Image::Pointer data, float recordTime,
   float BPHighPass, float BPLowPass,
   float alphaHighPass, float alphaLowPass)
 {
   bool powerOfTwo = false;
   int finalPower = 0;
   for (int i = 1; pow(2, i) <= data->GetDimension(1); ++i)
   {
     finalPower = i;
     if (pow(2, i) == data->GetDimension(1))
     {
       powerOfTwo = true;
     }
   }
   if (!powerOfTwo)
   {
     unsigned int dim[2] = { data->GetDimension(0), (unsigned int)pow(2,finalPower + 1) };
     data = ApplyResampling(data, dim);
   }
 
   MITK_INFO << data->GetDimension(0);
 
   // do a fourier transform, multiply with an appropriate window for the filter, and transform back
   typedef float PixelType;
   typedef itk::Image< PixelType, 3 > RealImageType;
   RealImageType::Pointer image;
 
   mitk::CastToItkImage(data, image);
 
   typedef itk::FFT1DRealToComplexConjugateImageFilter<RealImageType> ForwardFFTFilterType;
   typedef ForwardFFTFilterType::OutputImageType ComplexImageType;
   ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New();
   forwardFFTFilter->SetInput(image);
   forwardFFTFilter->SetDirection(1);
   try
   {
     forwardFFTFilter->UpdateOutputInformation();
   }
   catch (itk::ExceptionObject & error)
   {
     std::cerr << "Error: " << error << std::endl;
     MITK_WARN << "Bandpass could not be applied";
     return data;
   }
 
   float singleVoxel = 1 / (recordTime / data->GetDimension(1)) / 2 / 1000;
   float cutoffPixelHighPass = std::min(BPHighPass / singleVoxel, (float)data->GetDimension(1) / 2);
   float cutoffPixelLowPass = std::min(BPLowPass / singleVoxel, (float)data->GetDimension(1) / 2 - cutoffPixelHighPass);
 
   RealImageType::Pointer fftMultiplicator = BPFunction(data, cutoffPixelHighPass, cutoffPixelLowPass, alphaHighPass, alphaLowPass);
 
   typedef itk::MultiplyImageFilter< ComplexImageType,
     RealImageType,
     ComplexImageType >
     MultiplyFilterType;
   MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
   multiplyFilter->SetInput1(forwardFFTFilter->GetOutput());
   multiplyFilter->SetInput2(fftMultiplicator);
 
   /*itk::ComplexToModulusImageFilter<ComplexImageType, RealImageType>::Pointer toReal = itk::ComplexToModulusImageFilter<ComplexImageType, RealImageType>::New();
   toReal->SetInput(forwardFFTFilter->GetOutput());
   return GrabItkImageMemory(toReal->GetOutput());
   return GrabItkImageMemory(fftMultiplicator);   *///DEBUG
 
   typedef itk::FFT1DComplexConjugateToRealImageFilter< ComplexImageType, RealImageType > InverseFilterType;
   InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New();
   inverseFFTFilter->SetInput(multiplyFilter->GetOutput());
   inverseFFTFilter->SetDirection(1);
 
   return GrabItkImageMemory(inverseFFTFilter->GetOutput());
 }
 
 itk::Image<float, 3U>::Pointer mitk::PhotoacousticFilterService::BPFunction(mitk::Image::Pointer reference,
   int cutoffFrequencyPixelHighPass,
   int cutoffFrequencyPixelLowPass,
   float alphaHighPass, float alphaLowPass)
 {
   float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)];
 
   float width = reference->GetDimension(1) / 2.0 - (float)cutoffFrequencyPixelHighPass - (float)cutoffFrequencyPixelLowPass;
   float center = (float)cutoffFrequencyPixelHighPass / 2.0 + width / 2.0;
 
   for (unsigned int n = 0; n < reference->GetDimension(1); ++n)
   {
     imageData[reference->GetDimension(0)*n] = 0;
   }
   for (int n = 0; n < width; ++n)
   {
     imageData[reference->GetDimension(0)*n] = 1;
     if (n <= (alphaHighPass*(width - 1)) / 2.0)
     {
       if (alphaHighPass > 0.00001)
       {
         imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] =
           (1 + cos(itk::Math::pi*(2 * n / (alphaHighPass*(width - 1)) - 1))) / 2;
       }
       else
       {
         imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1;
       }
     }
     else if (n >= (width - 1)*(1 - alphaLowPass / 2)) //???
     {
       if (alphaLowPass > 0.00001)
       {
         imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] =
           (1 + cos(itk::Math::pi*(2 * n / (alphaLowPass*(width - 1)) + 1 - 2 / alphaLowPass))) / 2;
       }
       else
       {
         imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1;
       }
     }
     //MITK_INFO << "n:" << n << " is " << imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))];
   }
   MITK_INFO << "width: " << width << ", center: " << center << ", alphaHighPass: " << alphaHighPass << ", alphaLowPass: " << alphaLowPass;
 
   // mirror the first half of the image
   for (unsigned int n = reference->GetDimension(1) / 2; n < reference->GetDimension(1); ++n)
   {
     imageData[reference->GetDimension(0)*n] = imageData[(reference->GetDimension(1) - (n + 1)) * reference->GetDimension(0)];
   }
 
   // copy and paste to all lines
   for (unsigned int line = 1; line < reference->GetDimension(0); ++line)
   {
     for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample)
     {
       imageData[reference->GetDimension(0)*sample + line] = imageData[reference->GetDimension(0)*sample];
     }
   }
 
   typedef itk::Image< float, 3U >  ImageType;
   ImageType::RegionType region;
   ImageType::IndexType start;
   start.Fill(0);
 
   region.SetIndex(start);
 
   ImageType::SizeType size;
   size[0] = reference->GetDimension(0);
   size[1] = reference->GetDimension(1);
   size[2] = reference->GetDimension(2);
 
   region.SetSize(size);
 
   ImageType::SpacingType SpacingItk;
   SpacingItk[0] = reference->GetGeometry()->GetSpacing()[0];
   SpacingItk[1] = reference->GetGeometry()->GetSpacing()[1];
   SpacingItk[2] = reference->GetGeometry()->GetSpacing()[2];
 
   ImageType::Pointer image = ImageType::New();
   image->SetRegions(region);
   image->Allocate();
   image->FillBuffer(itk::NumericTraits<float>::Zero);
   image->SetSpacing(SpacingItk);
 
   ImageType::IndexType pixelIndex;
 
   for (unsigned int slice = 0; slice < reference->GetDimension(2); ++slice)
   {
     for (unsigned int line = 0; line < reference->GetDimension(0); ++line)
     {
       for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample)
       {
         pixelIndex[0] = line;
         pixelIndex[1] = sample;
         pixelIndex[2] = slice;
 
         image->SetPixel(pixelIndex, imageData[line + sample*reference->GetDimension(0)]);
       }
     }
   }
 
   delete[] imageData;
 
   return image;
 }