diff --git a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt index 7b8587334b..cc637416ca 100644 --- a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt +++ b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt @@ -1,16 +1,18 @@ set(dependencies_list MitkCore MitkAlgorithmsExt) IF(MITK_USE_OpenCL) add_definitions(-DPHOTOACOUSTICS_USE_GPU) set(dependencies_list ${dependencies_list} MitkOpenCL) message("Using OpenCL in PhotoacousticAlgorithms") ENDIF(MITK_USE_OpenCL) MITK_CREATE_MODULE( SUBPROJECTS DEPENDS ${dependencies_list} #AUTOLOAD_WITH MitkCore INCLUDE_DIRS PUBLIC Algorithms/ITKUltrasound Algorithms Algorithms/OCL INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} PACKAGE_DEPENDS ITK|ITKFFT+ITKImageCompose+ITKImageIntensity -) \ No newline at end of file +) + +add_subdirectory(test) diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index b3fdaff026..1cfe4eee02 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,25 +1,26 @@ file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*") set(CPP_FILES - source/mitkPhotoacousticImage.cpp + source/mitkPhotoacousticFilterService.cpp source/mitkPhotoacousticBeamformingFilter.cpp + source/mitkPhotoacousticBeamformingUtils.cpp source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp ) IF(MITK_USE_OpenCL) list(APPEND CPP_FILES source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp ) ENDIF(MITK_USE_OpenCL) set(RESOURCE_FILES BModeAbs.cl BModeAbsLog.cl UsedLinesCalculation.cl DelayCalculation.cl DMAS.cl DAS.cl sDMAS.cl ) diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h index 7908fccbae..84fc3a2bd2 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h @@ -1,151 +1,131 @@ /*=================================================================== 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. ===================================================================*/ #ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #define MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #include "mitkImageToImageFilter.h" #include #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "mitkPhotoacousticBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing an mitk::ImageToImageFilter for beamforming on both CPU and GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::BeamformingFilter::Configure(BeamformingSettings settings) * Whether the GPU is used can be set in the configuration. * For significant problems or important messages a string is written, which can be accessed via GetMessageString(). */ class BeamformingFilter : public ImageToImageFilter { public: mitkClassMacro(BeamformingFilter, ImageToImageFilter); itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** \brief Sets a new configuration to use * * @param settings The configuration set to use for beamforming */ void Configure(BeamformingSettings settings) { m_ConfOld = m_Conf; m_Conf = settings; } /** \brief Sets a new configuration to use * * The Filter writes important messages that can be retrieved through this method; if nothing is to be reported, it returns "noMessage". * @return The message */ std::string GetMessageString() { return m_Message; } /** \brief Sets a callback for progress checking * * An std::function can be set, through which progress of the currently updating filter is reported. * The integer argument is a number between 0 an 100 to indicate how far completion has been achieved, the std::string argument indicates what the filter is currently doing. */ void SetProgressHandle(std::function progressHandle); protected: BeamformingFilter(); ~BeamformingFilter() override; void GenerateInputRequestedRegion() override; void GenerateOutputInformation() override; void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; /** \brief The std::function, through which progress of the currently updating filter is reported. */ std::function m_ProgressHandle; /** \brief Pointer holding the Von-Hann apodization window for beamforming * @param samples the resolution at which the window is created */ float* VonHannFunction(int samples); /** \brief Function to create a Hamming apodization window * @param samples the resolution at which the window is created */ float* HammFunction(int samples); /** \brief Function to create a Box apodization window * @param samples the resolution at which the window is created */ float* BoxFunction(int samples); - - /** \brief Function to perform beamforming on CPU for a single line, using DAS and quadratic delay - */ - void DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); - /** \brief Function to perform beamforming on CPU for a single line, using DAS and spherical delay - */ - void DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); - /** \brief Function to perform beamforming on CPU for a single line, using DMAS and quadratic delay - */ - void DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); - /** \brief Function to perform beamforming on CPU for a single line, using DMAS and spherical delay - */ - void DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); - /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and quadratic delay - */ - void sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); - /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and spherical delay - */ - void sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); - float* m_OutputData; float* m_InputData; float* m_InputDataPuffer; /** \brief Pointer holding the Von-Hann apodization window for beamforming */ float* m_VonHannFunction; /** \brief Pointer holding the Hamming apodization window for beamforming */ float* m_HammFunction; /** \brief Pointer holding the Box apodization window for beamforming */ float* m_BoxFunction; /** \brief Current configuration set */ BeamformingSettings m_Conf; /** \brief Previous configuration set to selectively update used data for beamforming */ BeamformingSettings m_ConfOld; /** \brief Pointer to the GPU beamforming filter class; for performance reasons the filter is initialized within the constructor and kept for all later computations. */ mitk::PhotoacousticOCLBeamformingFilter::Pointer m_BeamformingOclFilter; /** \brief The message returned by mitk::BeamformingFilter::GetMessageString() */ std::string m_Message; }; } // namespace mitk #endif //MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h index a48b5e78a3..19d7f949c2 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h @@ -1,137 +1,137 @@ /*=================================================================== 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. ===================================================================*/ #ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_SETTINGS #define MITK_PHOTOACOUSTICS_BEAMFORMING_SETTINGS namespace mitk { /*! * \brief Class holding the configuration data for the beamforming filters mitk::BeamformingFilter and mitk::PhotoacousticOCLBeamformingFilter * * A detailed description can be seen below. All parameters should be set manually for successfull beamforming. */ class BeamformingSettings { public: /** \brief Pitch of the used transducer in [m]. */ float Pitch = 0.0003; /** \brief Speed of sound in the used medium in [m/s]. */ float SpeedOfSound = 1540; /** \brief This parameter is not neccessary to be set, as it's never used. */ float RecordTime = 0.00006; // [s] /** \brief The time spacing of the input image */ float TimeSpacing = 0.0000000000001; // [s] /** \brief The angle of the transducer elements */ float Angle = -1; /** \brief Flag whether processed image is a photoacoustic image or an ultrasound image */ bool isPhotoacousticImage = true; /** \brief How many transducer elements the used transducer had. */ unsigned short TransducerElements = 128; /** \brief How many vertical samples should be used in the final image. */ unsigned int SamplesPerLine = 2048; /** \brief How many lines should be reconstructed in the final image. */ unsigned int ReconstructionLines = 128; /** \brief Sets how many voxels should be cut off from the top of the image before beamforming, to potentially avoid artifacts. */ unsigned int upperCutoff = 0; /** \brief Sets whether only the slices selected by mitk::BeamformingSettings::CropBounds should be beamformed. */ bool partial = false; /** \brief Sets the first and last slice to be beamformed. */ unsigned int CropBounds[2] = { 0,0 }; /** \brief Sets the dimensions of the inputImage. */ unsigned int inputDim[3] = { 1,1,1 }; /** \brief Decides whether GPU computing should be used */ - bool UseGPU = true; + bool UseGPU = false; /** \brief Available delay calculation methods: * - Spherical delay for best results. * - A quadratic Taylor approximation for slightly faster results with hardly any quality loss. */ enum DelayCalc { QuadApprox, Spherical }; /** \brief Sets how the delays for beamforming should be calculated. */ DelayCalc DelayCalculationMethod = QuadApprox; /** \brief Available apodization functions: * - Hamming function. * - Von-Hann function. * - Box function. */ enum Apodization { Hamm, Hann, Box }; /** \brief Sets the used apodization function. */ Apodization Apod = Hann; /** \brief Sets the resolution of the apodization array (must be greater than 0). */ int apodizationArraySize = 128; /** \brief Available beamforming algorithms: * - DAS (Delay and sum). * - DMAS (Delay multiply and sum). */ enum BeamformingAlgorithm { DMAS, DAS, sDMAS}; /** \brief Sets the used beamforming algorithm. */ BeamformingAlgorithm Algorithm = DAS; /** \brief Sets whether after beamforming a bandpass should be automatically applied */ bool UseBP = false; /** \brief Sets the position at which lower frequencies are completely cut off in Hz. */ float BPHighPass = 50; /** \brief Sets the position at which higher frequencies are completely cut off in Hz. */ float BPLowPass = 50; /** \brief function for mitk::PhotoacousticOCLBeamformingFilter to check whether buffers need to be updated * this method only checks parameters relevant for the openCL implementation */ static bool SettingsChangedOpenCL(const BeamformingSettings& lhs, const BeamformingSettings& rhs) { return !((abs(lhs.Angle - rhs.Angle) < 0.01f) && // 0.01 degree error margin (lhs.Apod == rhs.Apod) && (lhs.DelayCalculationMethod == rhs.DelayCalculationMethod) && (lhs.isPhotoacousticImage == rhs.isPhotoacousticImage) && (abs(lhs.Pitch - rhs.Pitch) < 0.000001f) && // 0.0001 mm error margin (lhs.ReconstructionLines == rhs.ReconstructionLines) && (abs(lhs.RecordTime - rhs.RecordTime) < 0.00000001f) && // 10 ns error margin (lhs.SamplesPerLine == rhs.SamplesPerLine) && (abs(lhs.SpeedOfSound - rhs.SpeedOfSound) < 0.01f) && (abs(lhs.TimeSpacing - rhs.TimeSpacing) < 0.00000000001f) && //0.01 ns error margin (lhs.TransducerElements == rhs.TransducerElements)); } }; } #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingUtils.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingUtils.h new file mode 100644 index 0000000000..b243041b85 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingUtils.h @@ -0,0 +1,59 @@ +/*=================================================================== + +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. + +===================================================================*/ + +#ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER_UTILS +#define MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER_UTILS + +#include "mitkImageToImageFilter.h" +#include +#include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" +#include "mitkPhotoacousticBeamformingSettings.h" + +namespace mitk { + /*! + * \brief Class implementing util functionality for beamforming on CPU + * + */ + class PhotoacousticBeamformingUtils final + { + public: + /** \brief Function to perform beamforming on CPU for a single line, using DAS and quadratic delay + */ + static void DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config); + /** \brief Function to perform beamforming on CPU for a single line, using DAS and spherical delay + */ + static void DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config); + /** \brief Function to perform beamforming on CPU for a single line, using DMAS and quadratic delay + */ + static void DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config); + /** \brief Function to perform beamforming on CPU for a single line, using DMAS and spherical delay + */ + static void DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config); + /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and quadratic delay + */ + static void sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config); + /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and spherical delay + */ + static void sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config); + + protected: + PhotoacousticBeamformingUtils(); + + ~PhotoacousticBeamformingUtils(); + }; +} // namespace mitk + +#endif //MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER_UTILS diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h similarity index 93% rename from Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h rename to Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h index 504fdc1090..478b94ecd3 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h @@ -1,125 +1,125 @@ /*=================================================================== 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. ===================================================================*/ -#ifndef mitkPhotoacousticImage_H_HEADER_INCLUDED -#define mitkPhotoacousticImage_H_HEADER_INCLUDED +#ifndef mitkPhotoacousticFilterService_H_HEADER_INCLUDED +#define mitkPhotoacousticFilterService_H_HEADER_INCLUDED #include "itkObject.h" #include "mitkCommon.h" #include "mitkImage.h" #include #include "mitkPhotoacousticBeamformingSettings.h" #include "mitkPhotoacousticBeamformingFilter.h" #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { /*! * \brief Class holding methods to apply all Filters within the Photoacoustics Algorithms Module * * Implemented are: * - A B-Mode Filter * - A Resampling Filter * - Beamforming on GPU and CPU * - A Bandpass Filter */ - class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticImage : public itk::Object + class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticFilterService : public itk::Object { public: - mitkClassMacroItkParent(mitk::PhotoacousticImage, itk::Object); + mitkClassMacroItkParent(mitk::PhotoacousticFilterService, itk::Object); itkFactorylessNewMacro(Self); /** \brief Defines the methods for the B-Mode filter * Currently implemented are an Envelope Detection filter and a simple Absolute filter. */ enum BModeMethod { EnvelopeDetection, Abs }; /** \brief Applies a B-Mode Filter * * Applies a B-Mode filter using the given parameters. * @param inputImage The image to be processed. * @param method The kind of B-Mode Filter to be used. * @param UseGPU Setting this to true will allow the Filter to use the GPU. * @param UseLogFilter Setting this to true will apply a simple logarithm to the image after the B-Mode Filter has been applied. * @param resampleSpacing If this is set to 0, nothing will be done; otherwise, the image is resampled to a spacing of resampleSpacing mm per pixel. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseGPU = false, bool UseLogFilter = false, float resampleSpacing = 0.15); // mitk::Image::Pointer ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scatteringCoefficient); /** \brief Resamples the given image * * Resamples an image using the given parameters. * @param inputImage The image to be processed. * @param outputSize An array of dimensions the image should be resampled to. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]); /** \brief Beamforms the given image * * Resamples an image using the given parameters. * @param inputImage The image to be processed. * @param config The configuration set to be used for beamforming. * @param message A string into which potentially critical messages will be written. * @param progressHandle An std::function, through which progress of the currently updating filter is reported. * The integer argument is a number between 0 an 100 to indicate how far completion has been achieved, the std::string argument indicates what the filter is currently doing. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::string& message, std::function progressHandle = [](int, std::string) {}); /** \brief Crops the given image * * Crops an image in 3 dimension using the given parameters. * @param inputImage The image to be processed. * @param above How many voxels will be cut from the top of the image. * @param below How many voxels will be cut from the bottom of the image. * @param right How many voxels will be cut from the right side of the image. * @param left How many voxels will be cut from the left side of the image. * @param minSlice The first slice to be present in the resulting image. * @param maxSlice The last slice to be present in the resulting image. * @return The processed image is returned after the filter has finished. For the purposes of this module, the returned image is always of type float. */ mitk::Image::Pointer ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice); /** \brief Applies a Bandpass filter to the given image * * Applies a bandpass filter to the given image using the given parameters. * @param data The image to be processed. * @param recordTime The depth of the image in seconds. * @param BPHighPass The position at which Lower frequencies are completely cut off in Hz. * @param BPLowPass The position at which Higher frequencies are completely cut off in Hz. * @param alpha The tukey window parameter to control the shape of the bandpass filter: 0 will make it a Box function, 1 a Hann function. alpha can be set between those two bounds. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha); protected: - PhotoacousticImage(); - ~PhotoacousticImage() override; + PhotoacousticFilterService(); + ~PhotoacousticFilterService() override; /** \brief For performance reasons, an instance of the Beamforming filter is initialized as soon as possible and kept for all further uses. */ mitk::BeamformingFilter::Pointer m_BeamformingFilter; /** \brief Function that creates a Tukey function for the bandpass */ itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha); }; } // namespace mitk -#endif /* mitkPhotoacousticImage_H_HEADER_INCLUDED */ +#endif /* mitkPhotoacousticFilterService_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp deleted file mode 100644 index 343aecf11a..0000000000 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ /dev/null @@ -1,520 +0,0 @@ -/*=================================================================== - -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 "mitkPhotoacousticImage.h" -#include "ITKFilter/ITKUltrasound/itkBModeImageFilter.h" -#include "ITKFilter/itkPhotoacousticBModeImageFilter.h" -#include "mitkImageCast.h" -#include "mitkITKImageImport.h" -#include "mitkPhotoacousticBeamformingFilter.h" -#include -#include -#include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" - - -// itk dependencies -#include "itkImage.h" -#include "itkResampleImageFilter.h" -#include "itkCastImageFilter.h" -#include "itkCropImageFilter.h" -#include "itkRescaleIntensityImageFilter.h" -#include "itkIntensityWindowingImageFilter.h" -#include -#include "itkMultiplyImageFilter.h" -#include "itkBSplineInterpolateImageFunction.h" -#include - -// needed itk image filters -#include "mitkITKImageImport.h" -#include "itkFFTShiftImageFilter.h" -#include "itkMultiplyImageFilter.h" -#include "itkComplexToModulusImageFilter.h" -#include -#include "ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" -#include "ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" - -mitk::PhotoacousticImage::PhotoacousticImage() -{ - MITK_INFO << "[PhotoacousticImage Debug] created that image"; -} - -mitk::PhotoacousticImage::~PhotoacousticImage() -{ - MITK_INFO << "[PhotoacousticImage Debug] destroyed that image"; -} - -mitk::Image::Pointer mitk::PhotoacousticImage::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 TransformType; - - if (method == BModeMethod::Abs) - { - mitk::Image::Pointer input; - 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, inputImage->GetDimension(2) - 1); - - 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; - 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 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::ShapeDetection) - { - 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::PhotoacousticImage::ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scattering) -{ - typedef itk::Image< float, 3 > itkFloatImageType; - typedef itk::MultiplyImageFilter MultiplyImageFilterType; - - itkFloatImageType::Pointer itkImage; - mitk::CastToItkImage(inputImage, itkImage); - - MultiplyImageFilterType::Pointer multiplyFilter = MultiplyImageFilterType::New(); - multiplyFilter->SetInput1(itkImage); - multiplyFilter->SetInput2(m_FluenceCompResizedItk.at(m_ScatteringCoefficient)); - - return mitk::GrabItkImageMemory(multiplyFilter->GetOutput()); -}*/ - -mitk::Image::Pointer mitk::PhotoacousticImage::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 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(inputSizeItk[0]) / static_cast(outputSizeItk[0])); - outputSpacingItk[1] = itkImage->GetSpacing()[1] * (static_cast(inputSizeItk[1]) / static_cast(outputSizeItk[1])); - outputSpacingItk[2] = itkImage->GetSpacing()[2]; - - typedef itk::IdentityTransform 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::PhotoacousticImage::ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice) -{ - - unsigned int inputDim[3] = { inputImage->GetDimension(0), inputImage->GetDimension(1), inputImage->GetDimension(2) }; - unsigned int outputDim[3] = { inputImage->GetDimension(0) - left - right, inputImage->GetDimension(1) - (unsigned int)above - (unsigned int)below, (unsigned int)maxSlice - (unsigned int)minSlice + 1 }; - - void* inputData; - float* outputData = new float[outputDim[0] * outputDim[1] * outputDim[2]]; - - ImageReadAccessor acc(inputImage); - inputData = const_cast(acc.GetData()); - - // convert the data to float by default - // as of now only those float, short, float are used at all... though it's easy to add other ones - if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") - { - // copy the data into the cropped image - for (unsigned short sl = 0; sl < outputDim[2]; ++sl) - { - for (unsigned short l = 0; l < outputDim[0]; ++l) - { - for (unsigned short s = 0; s < outputDim[1]; ++s) - { - outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((float*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; - } - } - } - } - else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (short)" || inputImage->GetPixelType().GetTypeAsString() == " (short)") - { - // copy the data unsigned shorto the cropped image - for (unsigned short sl = 0; sl < outputDim[2]; ++sl) - { - for (unsigned short l = 0; l < outputDim[0]; ++l) - { - for (unsigned short s = 0; s < outputDim[1]; ++s) - { - outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((short*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; - } - } - } - } - else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (double)") - { - // copy the data unsigned shorto the cropped image - for (unsigned short sl = 0; sl < outputDim[2]; ++sl) - { - for (unsigned short l = 0; l < outputDim[0]; ++l) - { - for (unsigned short s = 0; s < outputDim[1]; ++s) - { - outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((double*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; - } - } - } - } - else - { - MITK_INFO << "Could not determine pixel type"; - } - - mitk::Image::Pointer output = mitk::Image::New(); - output->Initialize(mitk::MakeScalarPixelType(), 3, outputDim); - output->SetSpacing(inputImage->GetGeometry()->GetSpacing()); - output->SetImportVolume(outputData, 0, 0, mitk::Image::ReferenceMemory); - - return output; -} - -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, int cutoff, std::function progressHandle) -{ - config.RecordTime = config.RecordTime - (cutoff) / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping - progressHandle(0, "cropping image"); - if (!config.partial) - { - config.CropBounds[0] = 0; - config.CropBounds[1] = inputImage->GetDimension(2) - 1; - } - - Image::Pointer processedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); - config.inputDim[0] = processedImage->GetDimension(0); - config.inputDim[1] = processedImage->GetDimension(1); - config.inputDim[2] = processedImage->GetDimension(2); - - // perform the beamforming - BeamformingFilter::Pointer Beamformer = BeamformingFilter::New(); - Beamformer->SetInput(processedImage); - Beamformer->Configure(config); - Beamformer->SetProgressHandle(progressHandle); - Beamformer->UpdateLargestPossibleRegion(); - - processedImage = Beamformer->GetOutput(); - - return processedImage; -} - -mitk::Image::Pointer mitk::PhotoacousticImage::BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha) -{ - 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 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, alpha); - - typedef itk::MultiplyImageFilter< ComplexImageType, - RealImageType, - ComplexImageType > - MultiplyFilterType; - MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); - multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); - multiplyFilter->SetInput2(fftMultiplicator); - - /*itk::ComplexToModulusImageFilter::Pointer toReal = itk::ComplexToModulusImageFilter::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::Pointer mitk::PhotoacousticImage::BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha) -{ - float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; - - // tukey window - float width = reference->GetDimension(1) / 2 - (float)cutoffFrequencyPixelHighPass - (float)cutoffFrequencyPixelLowPass; - float center = (float)cutoffFrequencyPixelHighPass / 2 + width / 2; - - - MITK_INFO << width << "width " << center << "center " << alpha; - - for (unsigned int n = 0; n < reference->GetDimension(1); ++n) - { - imageData[reference->GetDimension(0)*n] = 0; - } - - for (int n = 0; n < width; ++n) - { - if (n <= (alpha*(width - 1)) / 2) - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(itk::Math::pi*(2 * n / (alpha*(width - 1)) - 1))) / 2; - } - else if (n >= (width - 1)*(1 - alpha / 2) && n <= (width - 1)) - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(itk::Math::pi*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; - } - else - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; - } - } - - // Butterworth-Filter - /* - // first, write the HighPass - if (cutoffFrequencyPixelHighPass != reference->GetDimension(1) / 2) - { - for (int n = 0; n < reference->GetDimension(1) / 2; ++n) - { - imageData[reference->GetDimension(0)*n] = 1 / (1 + pow( - (float)n / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelHighPass) - , 2 * butterworthOrder)); - } - } - else - { - for (int n = 0; n < reference->GetDimension(1) / 2; ++n) - { - imageData[reference->GetDimension(0)*n] = 1; - } - } - - // now, the LowPass - for (int n = 0; n < reference->GetDimension(1) / 2; ++n) - { - imageData[reference->GetDimension(0)*n] *= 1 / (1 + pow( - (float)(reference->GetDimension(1) / 2 - 1 - n) / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelLowPass) - , 2 * butterworthOrder)); - } - */ - - // 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::Zero); - image->SetSpacing(SpacingItk); - - ImageType::IndexType pixelIndex; - - for (ImageType::IndexValueType slice = 0; slice < reference->GetDimension(2); ++slice) - { - for (ImageType::IndexValueType line = 0; line < reference->GetDimension(0); ++line) - { - for (ImageType::IndexValueType 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; -} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h deleted file mode 100644 index 2200124e8d..0000000000 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h +++ /dev/null @@ -1,50 +0,0 @@ -/*=================================================================== - -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. - -===================================================================*/ - -#ifndef mitkPhotoacousticImage_H_HEADER_INCLUDED -#define mitkPhotoacousticImage_H_HEADER_INCLUDED - -#include "itkObject.h" -#include "mitkCommon.h" -#include "mitkImage.h" -#include - -#include "mitkPhotoacousticBeamformingFilter.h" -#include "MitkPhotoacousticsAlgorithmsExports.h" - -namespace mitk { - - class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticImage : public itk::Object - { - public: - mitkClassMacroItkParent(mitk::PhotoacousticImage, itk::Object); - itkFactorylessNewMacro(Self); - enum BModeMethod { ShapeDetection, Abs }; - mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseGPU = false, bool UseLogFilter = false, float resampleSpacing = 0.15); -// mitk::Image::Pointer ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scatteringCoefficient); - mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]); - mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); - mitk::Image::Pointer ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice); - mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha); - protected: - PhotoacousticImage(); - virtual ~PhotoacousticImage(); - - itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha); - }; -} // namespace mitk - -#endif /* mitkPhotoacousticImage_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp index 20abc2c982..32edd9b4d6 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp @@ -1,801 +1,367 @@ /*=================================================================== mitkPhotoacousticBeamformingFilter 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 "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include "mitkImageCast.h" #include "mitkPhotoacousticBeamformingFilter.h" +#include "mitkPhotoacousticBeamformingUtils.h" mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr), m_Message("noMessage") { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); m_VonHannFunction = VonHannFunction(m_Conf.apodizationArraySize); m_HammFunction = HammFunction(m_Conf.apodizationArraySize); m_BoxFunction = BoxFunction(m_Conf.apodizationArraySize); } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { delete[] m_VonHannFunction; delete[] m_HammFunction; delete[] m_BoxFunction; } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); //GenerateTimeInInputRegion(output, input); } void mitk::BeamformingFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; itkDebugMacro(<< "GenerateOutputInformation()"); unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf.Pitch * m_Conf.TransducerElements * 1000 / m_Conf.ReconstructionLines; spacing[1] = (m_Conf.TimeSpacing * m_Conf.inputDim[1]) / 2 * m_Conf.SpeedOfSound * 1000 / m_Conf.SamplesPerLine; spacing[2] = 1; output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::BeamformingFilter::GenerateData() { GenerateOutputInformation(); mitk::Image::Pointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; float* ApodWindow; if (m_ConfOld.apodizationArraySize != m_Conf.apodizationArraySize) { delete[] m_VonHannFunction; delete[] m_HammFunction; delete[] m_BoxFunction; m_VonHannFunction = VonHannFunction(m_Conf.apodizationArraySize); m_HammFunction = HammFunction(m_Conf.apodizationArraySize); m_BoxFunction = BoxFunction(m_Conf.apodizationArraySize); m_ConfOld = m_Conf; } // set the appropiate apodization window switch (m_Conf.Apod) { case BeamformingSettings::Apodization::Hann: ApodWindow = m_VonHannFunction; break; case BeamformingSettings::Apodization::Hamm: ApodWindow = m_HammFunction; break; case BeamformingSettings::Apodization::Box: ApodWindow = m_BoxFunction; break; default: ApodWindow = m_BoxFunction; break; } auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf.UseGPU) { int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); // first, we check whether the dara is float, other formats are unsupported if (input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)") { m_InputData = (float*)inputReadAccessor.GetData(); } else { MITK_INFO << "Pixel type is not float, abort"; return; } m_OutputData = new float[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; // fill the image with zeros for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { - threads[line] = std::thread(&BeamformingFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + threads[line] = std::thread(&PhotoacousticBeamformingUtils::DASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize, m_Conf); } } else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { - threads[line] = std::thread(&BeamformingFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + threads[line] = std::thread(&PhotoacousticBeamformingUtils::DASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize, m_Conf); } } } else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { - threads[line] = std::thread(&BeamformingFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + threads[line] = std::thread(&PhotoacousticBeamformingUtils::DMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize, m_Conf); } } else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { - threads[line] = std::thread(&BeamformingFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + threads[line] = std::thread(&PhotoacousticBeamformingUtils::DMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize, m_Conf); } } } else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::sDMAS) { if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { - threads[line] = std::thread(&BeamformingFilter::sDMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + threads[line] = std::thread(&PhotoacousticBeamformingUtils::sDMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize, m_Conf); } } else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { - threads[line] = std::thread(&BeamformingFilter::sDMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); + threads[line] = std::thread(&PhotoacousticBeamformingUtils::sDMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize, m_Conf); } } } // wait for all lines to finish for (short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } output->SetSlice(m_OutputData, i); if (i % progInterval == 0) m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; m_OutputData = nullptr; m_InputData = nullptr; } } #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN else { try { // first, we check whether the data is float, other formats are unsupported if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type is not float, abort"; return; } unsigned long availableMemory = m_BeamformingOclFilter->GetDeviceMemory(); unsigned int batchSize = 16; unsigned int batches = (unsigned int)((float)input->GetDimension(2)/batchSize) + (input->GetDimension(2)%batchSize > 0); unsigned int batchDim[] = { input->GetDimension(0), input->GetDimension(1), batchSize }; unsigned int batchDimLast[] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % batchSize }; // the following safeguard is probably only needed for absurdly small GPU memory for (batchSize = 16; (unsigned long)batchSize * ((unsigned long)(batchDim[0] * batchDim[1]) * 4 + // single input image (float) (unsigned long)(m_Conf.ReconstructionLines * m_Conf.SamplesPerLine) * 4) // single output image (float) > availableMemory - (unsigned long)(m_Conf.ReconstructionLines / 2 * m_Conf.SamplesPerLine) * 2 - // Delays buffer (unsigned short) (unsigned long)(m_Conf.ReconstructionLines * m_Conf.SamplesPerLine) * 3 * 2 - // UsedLines buffer (unsigned short) 50 * 1024 * 1024; // 50 MB buffer for local data, system purposes etc --batchSize) {} if (batchSize < 1) { MITK_ERROR << "device memory too small for GPU beamforming"; return; } mitk::ImageReadAccessor copy(input); for(unsigned int i = 0; i < batches; ++i) { m_ProgressHandle(input->GetDimension(2)/batches * i, "performing reconstruction"); mitk::Image::Pointer inputBatch = mitk::Image::New(); if(i == batches - 1 && (input->GetDimension(2)%batchSize > 0)) { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDimLast); m_Conf.inputDim[2] = batchDimLast[2]; } else { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDim); m_Conf.inputDim[2] = batchDim[2]; } inputBatch->SetSpacing(input->GetGeometry()->GetSpacing()); inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i])); m_BeamformingOclFilter->SetApodisation(ApodWindow, m_Conf.apodizationArraySize); m_BeamformingOclFilter->SetConfig(m_Conf); m_BeamformingOclFilter->SetInput(inputBatch); m_BeamformingOclFilter->Update(); void* out = m_BeamformingOclFilter->GetOutput(); for(unsigned int slice = 0; slice < m_Conf.inputDim[2]; ++slice) { output->SetImportSlice( &(((float*)out)[m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * slice]), batchSize * i + slice, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; float* dummyData = new float[m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * m_Conf.inputDim[2]]; output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); m_Message = "An openCL error occurred; all GPU operations in this and the next session may be corrupted."; } } #endif m_TimeOfHeaderInitialization.Modified(); auto end = std::chrono::high_resolution_clock::now(); MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } float* mitk::BeamformingFilter::VonHannFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = (1 - cos(2 * itk::Math::pi * n / (samples - 1))) / 2; } return ApodWindow; } float* mitk::BeamformingFilter::HammFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 0.54 - 0.46*cos(2 * itk::Math::pi*n / (samples - 1)); } return ApodWindow; } float* mitk::BeamformingFilter::BoxFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 1; } return ApodWindow; } - -void mitk::BeamformingFilter::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) -{ - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; - - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; - - short AddSample = 0; - short maxLine = 0; - short minLine = 0; - float delayMultiplicator = 0; - float l_i = 0; - float s_i = 0; - - float part = 0.07 * inputL; - float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); - float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; - float apod_mult = 1; - - short usedLines = (maxLine - minLine); - - //quadratic delay - l_i = line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (float)sample / outputS * inputS / 2; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); - - apod_mult = (float)apodArraySize / (float)usedLines; - - delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; - - for (short l_s = minLine; l_s < maxLine; ++l_s) - { - AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - m_Conf.isPhotoacousticImage)*s_i; - if (AddSample < inputS && AddSample >= 0) - output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; - else - --usedLines; - } - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; - } -} - -void mitk::BeamformingFilter::DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) -{ - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; - - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; - - short AddSample = 0; - short maxLine = 0; - short minLine = 0; - float l_i = 0; - float s_i = 0; - - float part = 0.07 * inputL; - float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); - float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; - float apod_mult = 1; - - short usedLines = (maxLine - minLine); - - //exact delay - - l_i = (float)line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (float)sample / outputS * inputS / 2; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); - - apod_mult = (float)apodArraySize / (float)usedLines; - - for (short l_s = minLine; l_s < maxLine; ++l_s) - { - AddSample = (int)sqrt( - pow(s_i, 2) - + - pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) - ) + (1 - m_Conf.isPhotoacousticImage)*s_i; - if (AddSample < inputS && AddSample >= 0) - output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; - else - --usedLines; - } - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; - } -} - - -void mitk::BeamformingFilter::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) -{ - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; - - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; - - short maxLine = 0; - short minLine = 0; - float delayMultiplicator = 0; - float l_i = 0; - float s_i = 0; - - float part = 0.07 * inputL; - float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); - float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; - float apod_mult = 1; - - float mult = 0; - short usedLines = (maxLine - minLine); - - //quadratic delay - l_i = line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = sample / outputS * inputS / 2; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); - - apod_mult = (float)apodArraySize / (float)usedLines; - - delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; - - //calculate the AddSamples beforehand to save some time - short* AddSample = new short[maxLine - minLine]; - for (short l_s = 0; l_s < maxLine - minLine; ++l_s) - { - AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*s_i; - } - - float s_1 = 0; - float s_2 = 0; - - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) - { - if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) - { - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) - { - if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) - { - s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; - s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; - output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); - } - } - } - else - --usedLines; - } - - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); - - delete[] AddSample; - } -} - -void mitk::BeamformingFilter::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) -{ - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; - - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; - - short maxLine = 0; - short minLine = 0; - float l_i = 0; - float s_i = 0; - - float part = 0.07 * inputL; - float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); - float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; - float apod_mult = 1; - - float mult = 0; - - short usedLines = (maxLine - minLine); - - //exact delay - - l_i = (float)line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (float)sample / outputS * inputS / 2; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); - - apod_mult = (float)apodArraySize / (float)usedLines; - - //calculate the AddSamples beforehand to save some time - short* AddSample = new short[maxLine - minLine]; - for (short l_s = 0; l_s < maxLine - minLine; ++l_s) - { - AddSample[l_s] = (short)sqrt( - pow(s_i, 2) - + - pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) - ) + (1 - m_Conf.isPhotoacousticImage)*s_i; - } - - float s_1 = 0; - float s_2 = 0; - - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) - { - if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) - { - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) - { - if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) - { - s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; - s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; - output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); - } - } - } - else - --usedLines; - } - - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); - - delete[] AddSample; - } -} - -void mitk::BeamformingFilter::sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) -{ - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; - - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; - - short maxLine = 0; - short minLine = 0; - float delayMultiplicator = 0; - float l_i = 0; - float s_i = 0; - - float part = 0.07 * inputL; - float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); - float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; - float apod_mult = 1; - - float mult = 0; - short usedLines = (maxLine - minLine); - - //quadratic delay - l_i = line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = sample / outputS * inputS / 2; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); - - apod_mult = (float)apodArraySize / (float)usedLines; - - delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; - - //calculate the AddSamples beforehand to save some time - short* AddSample = new short[maxLine - minLine]; - for (short l_s = 0; l_s < maxLine - minLine; ++l_s) - { - AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*s_i; - } - - float s_1 = 0; - float s_2 = 0; - float sign = 0; - - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) - { - if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) - { - s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - sign += s_1; - - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) - { - if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) - { - s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; - - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; - output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); - } - } - } - else - --usedLines; - } - - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); - - delete[] AddSample; - } -} - -void mitk::BeamformingFilter::sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) -{ - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; - - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; - - short maxLine = 0; - short minLine = 0; - float l_i = 0; - float s_i = 0; - - float part = 0.07 * inputL; - float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); - float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; - float apod_mult = 1; - - float mult = 0; - - short usedLines = (maxLine - minLine); - - //exact delay - - l_i = (float)line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (float)sample / outputS * inputS / 2; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); - - apod_mult = (float)apodArraySize / (float)usedLines; - - //calculate the AddSamples beforehand to save some time - short* AddSample = new short[maxLine - minLine]; - for (short l_s = 0; l_s < maxLine - minLine; ++l_s) - { - AddSample[l_s] = (short)sqrt( - pow(s_i, 2) - + - pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) - ) + (1 - m_Conf.isPhotoacousticImage)*s_i; - } - - float s_1 = 0; - float s_2 = 0; - float sign = 0; - - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) - { - if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) - { - s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - sign += s_1; - - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) - { - if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) - { - s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; - - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; - output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); - } - } - } - else - --usedLines; - } - - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); - - delete[] AddSample; - } -} diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingUtils.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingUtils.cpp new file mode 100644 index 0000000000..59eea2d0b6 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingUtils.cpp @@ -0,0 +1,468 @@ +/*=================================================================== +mitkPhotoacousticBeamformingFilter +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 "mitkProperties.h" +#include "mitkImageReadAccessor.h" +#include +#include +#include +#include +#include +#include "mitkImageCast.h" +#include "mitkPhotoacousticBeamformingUtils.h" + +mitk::PhotoacousticBeamformingUtils::PhotoacousticBeamformingUtils() +{ +} + +mitk::PhotoacousticBeamformingUtils::~PhotoacousticBeamformingUtils() +{ +} + +void mitk::PhotoacousticBeamformingUtils::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short AddSample = 0; + short maxLine = 0; + short minLine = 0; + float delayMultiplicator = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(config.Angle / 360 * 2 * itk::Math::pi); + float part_multiplicator = tan_phi * config.TimeSpacing * config.SpeedOfSound / config.Pitch * inputL / config.TransducerElements; + float apod_mult = 1; + + short usedLines = (maxLine - minLine); + + //quadratic delay + l_i = line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + delayMultiplicator = pow((1 / (config.TimeSpacing*config.SpeedOfSound) * (config.Pitch*config.TransducerElements) / inputL), 2) / s_i / 2; + + for (short l_s = minLine; l_s < maxLine; ++l_s) + { + AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - config.isPhotoacousticImage)*s_i; + if (AddSample < inputS && AddSample >= 0) + output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; + else + --usedLines; + } + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; + } +} + +void mitk::PhotoacousticBeamformingUtils::DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short AddSample = 0; + short maxLine = 0; + short minLine = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(config.Angle / 360 * 2 * itk::Math::pi); + float part_multiplicator = tan_phi * config.TimeSpacing * config.SpeedOfSound / config.Pitch * inputL / (float)config.TransducerElements; + float apod_mult = 1; + + short usedLines = (maxLine - minLine); + + //exact delay + + l_i = (float)line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + for (short l_s = minLine; l_s < maxLine; ++l_s) + { + AddSample = (int)sqrt( + pow(s_i, 2) + + + pow((1 / (config.TimeSpacing*config.SpeedOfSound) * (((float)l_s - l_i)*config.Pitch*(float)config.TransducerElements) / inputL), 2) + ) + (1 - config.isPhotoacousticImage)*s_i; + if (AddSample < inputS && AddSample >= 0) + output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; + else + --usedLines; + } + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; + } +} + + +void mitk::PhotoacousticBeamformingUtils::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short maxLine = 0; + short minLine = 0; + float delayMultiplicator = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(config.Angle / 360 * 2 * itk::Math::pi); + float part_multiplicator = tan_phi * config.TimeSpacing * config.SpeedOfSound / config.Pitch * inputL / (float)config.TransducerElements; + float apod_mult = 1; + + float mult = 0; + short usedLines = (maxLine - minLine); + + //quadratic delay + l_i = line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + delayMultiplicator = pow((1 / (config.TimeSpacing*config.SpeedOfSound) * (config.Pitch*config.TransducerElements) / inputL), 2) / s_i / 2; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) + { + AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - config.isPhotoacousticImage)*s_i; + } + + float s_1 = 0; + float s_2 = 0; + + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + { + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + { + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); + + delete[] AddSample; + } +} + +void mitk::PhotoacousticBeamformingUtils::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short maxLine = 0; + short minLine = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(config.Angle / 360 * 2 * itk::Math::pi); + float part_multiplicator = tan_phi * config.TimeSpacing * config.SpeedOfSound / config.Pitch * inputL / (float)config.TransducerElements; + float apod_mult = 1; + + float mult = 0; + + short usedLines = (maxLine - minLine); + + //exact delay + + l_i = (float)line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) + { + AddSample[l_s] = (short)sqrt( + pow(s_i, 2) + + + pow((1 / (config.TimeSpacing*config.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*config.Pitch*(float)config.TransducerElements) / inputL), 2) + ) + (1 - config.isPhotoacousticImage)*s_i; + } + + float s_1 = 0; + float s_2 = 0; + + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + { + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + { + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); + + delete[] AddSample; + } +} + +void mitk::PhotoacousticBeamformingUtils::sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short maxLine = 0; + short minLine = 0; + float delayMultiplicator = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(config.Angle / 360 * 2 * itk::Math::pi); + float part_multiplicator = tan_phi * config.TimeSpacing * config.SpeedOfSound / config.Pitch * inputL / (float)config.TransducerElements; + float apod_mult = 1; + + float mult = 0; + short usedLines = (maxLine - minLine); + + //quadratic delay + l_i = line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + delayMultiplicator = pow((1 / (config.TimeSpacing*config.SpeedOfSound) * (config.Pitch*config.TransducerElements) / inputL), 2) / s_i / 2; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) + { + AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - config.isPhotoacousticImage)*s_i; + } + + float s_1 = 0; + float s_2 = 0; + float sign = 0; + + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + { + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + { + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + sign += s_1; + + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); + + delete[] AddSample; + } +} + +void mitk::PhotoacousticBeamformingUtils::sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize, const mitk::BeamformingSettings& config) +{ + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short maxLine = 0; + short minLine = 0; + float l_i = 0; + float s_i = 0; + + float part = 0.07 * inputL; + float tan_phi = std::tan(config.Angle / 360 * 2 * itk::Math::pi); + float part_multiplicator = tan_phi * config.TimeSpacing * config.SpeedOfSound / config.Pitch * inputL / (float)config.TransducerElements; + float apod_mult = 1; + + float mult = 0; + + short usedLines = (maxLine - minLine); + + //exact delay + + l_i = (float)line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * inputS / 2; + + part = part_multiplicator*s_i; + + if (part < 1) + part = 1; + + maxLine = (short)std::min((l_i + part) + 1, inputL); + minLine = (short)std::max((l_i - part), 0.0f); + usedLines = (maxLine - minLine); + + apod_mult = (float)apodArraySize / (float)usedLines; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) + { + AddSample[l_s] = (short)sqrt( + pow(s_i, 2) + + + pow((1 / (config.TimeSpacing*config.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*config.Pitch*(float)config.TransducerElements) / inputL), 2) + ) + (1 - config.isPhotoacousticImage)*s_i; + } + + float s_1 = 0; + float s_2 = 0; + float sign = 0; + + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + { + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + { + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + sign += s_1; + + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); + + delete[] AddSample; + } +} diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticFilterService.cpp similarity index 92% rename from Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp rename to Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticFilterService.cpp index d4961e2260..e5b3eb8d51 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticFilterService.cpp @@ -1,539 +1,539 @@ /*=================================================================== 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 "mitkPhotoacousticImage.h" +#include "mitkPhotoacousticFilterService.h" #include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" #include "../ITKFilter/itkPhotoacousticBModeImageFilter.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" #include "mitkPhotoacousticBeamformingFilter.h" #include #include #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" // itk dependencies #include "itkImage.h" #include "itkResampleImageFilter.h" #include "itkCastImageFilter.h" #include "itkCropImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include "itkIntensityWindowingImageFilter.h" #include #include "itkMultiplyImageFilter.h" #include "itkBSplineInterpolateImageFunction.h" #include // needed itk image filters #include "mitkITKImageImport.h" #include "itkFFTShiftImageFilter.h" #include "itkMultiplyImageFilter.h" #include "itkComplexToModulusImageFilter.h" #include #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" -mitk::PhotoacousticImage::PhotoacousticImage() : m_BeamformingFilter(BeamformingFilter::New()) +mitk::PhotoacousticFilterService::PhotoacousticFilterService() : m_BeamformingFilter(BeamformingFilter::New()) { - MITK_INFO << "[PhotoacousticImage Debug] created that image"; + MITK_INFO << "[PhotoacousticFilterService Debug] created that image"; } -mitk::PhotoacousticImage::~PhotoacousticImage() +mitk::PhotoacousticFilterService::~PhotoacousticFilterService() { - MITK_INFO << "[PhotoacousticImage Debug] destroyed that image"; + MITK_INFO << "[PhotoacousticFilterService Debug] destroyed that image"; } -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseGPU, bool UseLogFilter, float resampleSpacing) +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 TransformType; if (method == BModeMethod::Abs) { mitk::Image::Pointer input; 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, inputImage->GetDimension(2) - 1); 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; 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 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::PhotoacousticImage::ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scattering) +/*mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scattering) { typedef itk::Image< float, 3 > itkFloatImageType; typedef itk::MultiplyImageFilter MultiplyImageFilterType; itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); MultiplyImageFilterType::Pointer multiplyFilter = MultiplyImageFilterType::New(); multiplyFilter->SetInput1(itkImage); multiplyFilter->SetInput2(m_FluenceCompResizedItk.at(m_ScatteringCoefficient)); return mitk::GrabItkImageMemory(multiplyFilter->GetOutput()); }*/ -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]) +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 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(inputSizeItk[0]) / static_cast(outputSizeItk[0])); outputSpacingItk[1] = itkImage->GetSpacing()[1] * (static_cast(inputSizeItk[1]) / static_cast(outputSizeItk[1])); outputSpacingItk[2] = itkImage->GetSpacing()[2]; typedef itk::IdentityTransform 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::PhotoacousticImage::ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice) { unsigned int inputDim[3] = { inputImage->GetDimension(0), inputImage->GetDimension(1), inputImage->GetDimension(2) }; unsigned int outputDim[3] = { inputImage->GetDimension(0) - left - right, inputImage->GetDimension(1) - (unsigned int)above - (unsigned int)below, (unsigned int)maxSlice - (unsigned int)minSlice + 1 }; void* inputData; float* outputData = new float[outputDim[0] * outputDim[1] * outputDim[2]]; ImageReadAccessor acc(inputImage); inputData = const_cast(acc.GetData()); // convert the data to float by default // as of now only float, short, double are used at all. if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") { // copy the data into the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((float*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (short)" || inputImage->GetPixelType().GetTypeAsString() == " (short)") { // copy the data to the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((short*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (double)") { // copy the data to the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((double*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else { MITK_INFO << "Could not determine pixel type"; } mitk::Image::Pointer output = mitk::Image::New(); output->Initialize(mitk::MakeScalarPixelType(), 3, outputDim); output->SetSpacing(inputImage->GetGeometry()->GetSpacing()); output->SetImportVolume(outputData, 0, 0, mitk::Image::ReferenceMemory); return output; } -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::string& message, std::function progressHandle) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::string& message, std::function progressHandle) { Image::Pointer processedImage = inputImage; if (inputImage->GetDimension() != 3) { processedImage->Initialize(mitk::MakeScalarPixelType(), 3, inputImage->GetDimensions()); processedImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); mitk::ImageReadAccessor copy(inputImage); processedImage->SetImportVolume(copy.GetData()); } config.RecordTime = config.RecordTime - (float)(config.upperCutoff) / (float)inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping progressHandle(0, "converting image"); if (!config.partial) { config.CropBounds[0] = 0; config.CropBounds[1] = inputImage->GetDimension(2) - 1; } processedImage = ApplyCropping(inputImage, config.upperCutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); config.inputDim[0] = processedImage->GetDimension(0); config.inputDim[1] = processedImage->GetDimension(1); config.inputDim[2] = processedImage->GetDimension(2); // perform the beamforming m_BeamformingFilter->SetInput(processedImage); m_BeamformingFilter->Configure(config); m_BeamformingFilter->SetProgressHandle(progressHandle); m_BeamformingFilter->UpdateLargestPossibleRegion(); processedImage = m_BeamformingFilter->GetOutput(); message = m_BeamformingFilter->GetMessageString(); return processedImage; } -mitk::Image::Pointer mitk::PhotoacousticImage::BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha) +mitk::Image::Pointer mitk::PhotoacousticFilterService::BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha) { 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 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, alpha); typedef itk::MultiplyImageFilter< ComplexImageType, RealImageType, ComplexImageType > MultiplyFilterType; MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); multiplyFilter->SetInput2(fftMultiplicator); /*itk::ComplexToModulusImageFilter::Pointer toReal = itk::ComplexToModulusImageFilter::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::Pointer mitk::PhotoacousticImage::BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha) +itk::Image::Pointer mitk::PhotoacousticFilterService::BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha) { float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; // tukey window float width = reference->GetDimension(1) / 2 - (float)cutoffFrequencyPixelHighPass - (float)cutoffFrequencyPixelLowPass; float center = (float)cutoffFrequencyPixelHighPass / 2 + width / 2; MITK_INFO << width << "width " << center << "center " << alpha; for (unsigned int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = 0; } if (alpha < 0.00001) { for (int n = 0; n < width; ++n) { if (n <= (alpha*(width - 1)) / 2) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(itk::Math::pi*(2 * n / (alpha*(width - 1)) - 1))) / 2; } else if (n >= (width - 1)*(1 - alpha / 2)) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(itk::Math::pi*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; } else { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; } } } else { for (int n = 0; n < width; ++n) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; } } // Butterworth-Filter /* // first, write the HighPass if (cutoffFrequencyPixelHighPass != reference->GetDimension(1) / 2) { for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] = 1 / (1 + pow( (float)n / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelHighPass) , 2 * butterworthOrder)); } } else { for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] = 1; } } // now, the LowPass for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] *= 1 / (1 + pow( (float)(reference->GetDimension(1) / 2 - 1 - n) / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelLowPass) , 2 * butterworthOrder)); } */ // 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::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; } diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkphotoacousticbeamformingfilterutils.cpp similarity index 50% copy from Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp copy to Modules/PhotoacousticsAlgorithms/source/mitkphotoacousticbeamformingfilterutils.cpp index 20abc2c982..b447706ece 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkphotoacousticbeamformingfilterutils.cpp @@ -1,801 +1,468 @@ /*=================================================================== mitkPhotoacousticBeamformingFilter 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 "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include "mitkImageCast.h" -#include "mitkPhotoacousticBeamformingFilter.h" +#include "mitkPhotoacousticBeamformingUtils.h" -mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr), m_Message("noMessage") +mitk::PhotoacousticBeamformingUtils::PhotoacousticBeamformingUtils() { - this->SetNumberOfIndexedInputs(1); - this->SetNumberOfRequiredInputs(1); - - m_ProgressHandle = [](int, std::string) {}; - - m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); - - m_VonHannFunction = VonHannFunction(m_Conf.apodizationArraySize); - m_HammFunction = HammFunction(m_Conf.apodizationArraySize); - m_BoxFunction = BoxFunction(m_Conf.apodizationArraySize); -} - -void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) -{ - m_ProgressHandle = progressHandle; -} - -mitk::BeamformingFilter::~BeamformingFilter() -{ - delete[] m_VonHannFunction; - delete[] m_HammFunction; - delete[] m_BoxFunction; -} - -void mitk::BeamformingFilter::GenerateInputRequestedRegion() -{ - Superclass::GenerateInputRequestedRegion(); - - mitk::Image* output = this->GetOutput(); - mitk::Image* input = const_cast (this->GetInput()); - if (!output->IsInitialized()) - { - return; - } - - input->SetRequestedRegionToLargestPossibleRegion(); - //GenerateTimeInInputRegion(output, input); -} - -void mitk::BeamformingFilter::GenerateOutputInformation() -{ - mitk::Image::ConstPointer input = this->GetInput(); - mitk::Image::Pointer output = this->GetOutput(); - - if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) - return; - - itkDebugMacro(<< "GenerateOutputInformation()"); - - unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) }; - output->Initialize(mitk::MakeScalarPixelType(), 3, dim); - - mitk::Vector3D spacing; - spacing[0] = m_Conf.Pitch * m_Conf.TransducerElements * 1000 / m_Conf.ReconstructionLines; - spacing[1] = (m_Conf.TimeSpacing * m_Conf.inputDim[1]) / 2 * m_Conf.SpeedOfSound * 1000 / m_Conf.SamplesPerLine; - spacing[2] = 1; - - output->GetGeometry()->SetSpacing(spacing); - output->GetGeometry()->Modified(); - output->SetPropertyList(input->GetPropertyList()->Clone()); - - m_TimeOfHeaderInitialization.Modified(); } - -void mitk::BeamformingFilter::GenerateData() +mitk::PhotoacousticBeamformingUtils::~PhotoacousticBeamformingUtils() { - GenerateOutputInformation(); - mitk::Image::Pointer input = this->GetInput(); - mitk::Image::Pointer output = this->GetOutput(); - - if (!output->IsInitialized()) - return; - - float* ApodWindow; - if (m_ConfOld.apodizationArraySize != m_Conf.apodizationArraySize) - { - delete[] m_VonHannFunction; - delete[] m_HammFunction; - delete[] m_BoxFunction; - m_VonHannFunction = VonHannFunction(m_Conf.apodizationArraySize); - m_HammFunction = HammFunction(m_Conf.apodizationArraySize); - m_BoxFunction = BoxFunction(m_Conf.apodizationArraySize); - - m_ConfOld = m_Conf; - } - - // set the appropiate apodization window - switch (m_Conf.Apod) - { - case BeamformingSettings::Apodization::Hann: - ApodWindow = m_VonHannFunction; - break; - case BeamformingSettings::Apodization::Hamm: - ApodWindow = m_HammFunction; - break; - case BeamformingSettings::Apodization::Box: - ApodWindow = m_BoxFunction; - break; - default: - ApodWindow = m_BoxFunction; - break; - } - - auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... - - if (!m_Conf.UseGPU) - { - int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; - // the interval at which we update the gui progress bar - - float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; - float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; - - for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied - { - mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); - - // first, we check whether the dara is float, other formats are unsupported - if (input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)") - { - m_InputData = (float*)inputReadAccessor.GetData(); - } - else - { - MITK_INFO << "Pixel type is not float, abort"; - return; - } - - m_OutputData = new float[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; - - // fill the image with zeros - for (int l = 0; l < outputDim[0]; ++l) - { - for (int s = 0; s < outputDim[1]; ++s) - { - m_OutputData[l*(short)outputDim[1] + s] = 0; - } - } - - std::thread *threads = new std::thread[(short)outputDim[0]]; - - // every line will be beamformed in a seperate thread - if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DAS) - { - if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) - { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); - } - } - else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) - { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); - } - } - } - else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DMAS) - { - if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) - { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); - } - } - else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) - { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); - } - } - } - else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::sDMAS) - { - if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::QuadApprox) - { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingFilter::sDMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); - } - } - else if (m_Conf.DelayCalculationMethod == BeamformingSettings::DelayCalc::Spherical) - { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingFilter::sDMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, m_Conf.apodizationArraySize); - } - } - } - // wait for all lines to finish - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line].join(); - } - - output->SetSlice(m_OutputData, i); - - if (i % progInterval == 0) - m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction"); - - delete[] m_OutputData; - m_OutputData = nullptr; - m_InputData = nullptr; - } - } - #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN - else - { - try - { - // first, we check whether the data is float, other formats are unsupported - if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) - { - MITK_ERROR << "Pixel type is not float, abort"; - return; - } - - unsigned long availableMemory = m_BeamformingOclFilter->GetDeviceMemory(); - - unsigned int batchSize = 16; - unsigned int batches = (unsigned int)((float)input->GetDimension(2)/batchSize) + (input->GetDimension(2)%batchSize > 0); - - unsigned int batchDim[] = { input->GetDimension(0), input->GetDimension(1), batchSize }; - unsigned int batchDimLast[] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % batchSize }; - - // the following safeguard is probably only needed for absurdly small GPU memory - for (batchSize = 16; - (unsigned long)batchSize * - ((unsigned long)(batchDim[0] * batchDim[1]) * 4 + // single input image (float) - (unsigned long)(m_Conf.ReconstructionLines * m_Conf.SamplesPerLine) * 4) // single output image (float) - > availableMemory - - (unsigned long)(m_Conf.ReconstructionLines / 2 * m_Conf.SamplesPerLine) * 2 - // Delays buffer (unsigned short) - (unsigned long)(m_Conf.ReconstructionLines * m_Conf.SamplesPerLine) * 3 * 2 - // UsedLines buffer (unsigned short) - 50 * 1024 * 1024; // 50 MB buffer for local data, system purposes etc - --batchSize) - {} - if (batchSize < 1) - { - MITK_ERROR << "device memory too small for GPU beamforming"; - return; - } - - mitk::ImageReadAccessor copy(input); - - for(unsigned int i = 0; i < batches; ++i) - { - m_ProgressHandle(input->GetDimension(2)/batches * i, "performing reconstruction"); - - mitk::Image::Pointer inputBatch = mitk::Image::New(); - if(i == batches - 1 && (input->GetDimension(2)%batchSize > 0)) - { - inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDimLast); - m_Conf.inputDim[2] = batchDimLast[2]; - } - else - { - inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDim); - m_Conf.inputDim[2] = batchDim[2]; - } - - inputBatch->SetSpacing(input->GetGeometry()->GetSpacing()); - - inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i])); - - m_BeamformingOclFilter->SetApodisation(ApodWindow, m_Conf.apodizationArraySize); - m_BeamformingOclFilter->SetConfig(m_Conf); - m_BeamformingOclFilter->SetInput(inputBatch); - m_BeamformingOclFilter->Update(); - - void* out = m_BeamformingOclFilter->GetOutput(); - - for(unsigned int slice = 0; slice < m_Conf.inputDim[2]; ++slice) - { - output->SetImportSlice( - &(((float*)out)[m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * slice]), - batchSize * i + slice, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); - } - } - } - catch (mitk::Exception &e) - { - std::string errorMessage = "Caught unexpected exception "; - errorMessage.append(e.what()); - MITK_ERROR << errorMessage; - - float* dummyData = new float[m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * m_Conf.inputDim[2]]; - output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); - - m_Message = "An openCL error occurred; all GPU operations in this and the next session may be corrupted."; - } - } - #endif - m_TimeOfHeaderInitialization.Modified(); - - auto end = std::chrono::high_resolution_clock::now(); - MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; -} - -float* mitk::BeamformingFilter::VonHannFunction(int samples) -{ - float* ApodWindow = new float[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = (1 - cos(2 * itk::Math::pi * n / (samples - 1))) / 2; - } - - return ApodWindow; -} - -float* mitk::BeamformingFilter::HammFunction(int samples) -{ - float* ApodWindow = new float[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = 0.54 - 0.46*cos(2 * itk::Math::pi*n / (samples - 1)); - } - - return ApodWindow; -} - -float* mitk::BeamformingFilter::BoxFunction(int samples) -{ - float* ApodWindow = new float[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = 1; - } - - return ApodWindow; } -void mitk::BeamformingFilter::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +void mitk::PhotoacousticBeamformingUtils::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - m_Conf.isPhotoacousticImage)*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } -void mitk::BeamformingFilter::DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +void mitk::PhotoacousticBeamformingUtils::DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = (int)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.isPhotoacousticImage)*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } -void mitk::BeamformingFilter::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +void mitk::PhotoacousticBeamformingUtils::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*s_i; } float s_1 = 0; float s_2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } -void mitk::BeamformingFilter::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +void mitk::PhotoacousticBeamformingUtils::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.isPhotoacousticImage)*s_i; } float s_1 = 0; float s_2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } -void mitk::BeamformingFilter::sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +void mitk::PhotoacousticBeamformingUtils::sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*s_i; } float s_1 = 0; float s_2 = 0; float sign = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; sign += s_1; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); delete[] AddSample; } } -void mitk::BeamformingFilter::sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) +void mitk::PhotoacousticBeamformingUtils::sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / (float)m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (((float)minLine + (float)l_s - l_i)*m_Conf.Pitch*(float)m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.isPhotoacousticImage)*s_i; } float s_1 = 0; float s_2 = 0; float sign = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; sign += s_1; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); delete[] AddSample; } } diff --git a/Modules/PhotoacousticsAlgorithms/test/CMakeLists.txt b/Modules/PhotoacousticsAlgorithms/test/CMakeLists.txt new file mode 100644 index 0000000000..153cd81e2e --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/test/CMakeLists.txt @@ -0,0 +1 @@ +MITK_CREATE_MODULE_TESTS() diff --git a/Modules/PhotoacousticsAlgorithms/test/files.cmake b/Modules/PhotoacousticsAlgorithms/test/files.cmake new file mode 100644 index 0000000000..6e42fdeff9 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/test/files.cmake @@ -0,0 +1,21 @@ +set(MODULE_TESTS + # IMPORTANT: If you plan to deactivate / comment out a test please write a bug number to the commented out line of code. + # + # Example: #mitkMyTest #this test is commented out because of bug 12345 + # + # It is important that the bug is open and that the test will be activated again before the bug is closed. This assures that + # no test is forgotten after it was commented out. If there is no bug for your current problem, please add a new one and + # mark it as critical. + + ################## ON THE FENCE TESTS ################################################# + # none + + ################## DISABLED TESTS ##################################################### + + ################# RUNNING TESTS ####################################################### + + + mitkPAFilterServiceTest.cpp + +) +set(RESOURCE_FILES) diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp new file mode 100644 index 0000000000..571018e933 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp @@ -0,0 +1,89 @@ +/*=================================================================== + +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 +#include + +#include +#include +#include +#include + +// us +#include +#include +#include +#include +#include + +#include +#include + +class mitkPAFilterServiceTestSuite : public mitk::TestFixture +{ + CPPUNIT_TEST_SUITE(mitkPAFilterServiceTestSuite); + MITK_TEST(testRunning); + CPPUNIT_TEST_SUITE_END(); + +private: + + mitk::PhotoacousticFilterService::Pointer m_PhotoacousticFilterService; + +public: + + void setUp() override + { + m_PhotoacousticFilterService = mitk::PhotoacousticFilterService::New(); + } + + void testRunning() + { + int xDim = 3; + int yDim = 3; + int length = yDim * xDim; + int* testArray = new int [length]; + for (int i = 0; i < length; ++i) + { + testArray[i] = 0; + } + unsigned int NUMBER_OF_SPATIAL_DIMENSIONS = 2; + auto* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; + dimensions[0] = yDim; + dimensions[1] = xDim; + mitk::PixelType pixelType = mitk::MakeScalarPixelType(); + mitk::Image::Pointer testImage = mitk::Image::New(); + testImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); + testImage->SetImportSlice(testArray, 0, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); + std::string message = "Es ist egal. Das taucht nicht auf!"; + mitk::BeamformingSettings config; + auto output = m_PhotoacousticFilterService->ApplyBeamforming(testImage, config, message); + + mitk::ImageReadAccessor readAccess(output, output->GetVolumeData()); + const int* outputArray = (const int*)readAccess.GetData(); + + for (int i = 0; i < length; ++i) + { + CPPUNIT_ASSERT(outputArray[i]==0); + } + } + + void tearDown() override + { + m_PhotoacousticFilterService = nullptr; + } +}; + +MITK_TEST_SUITE_REGISTRATION(mitkPAFilterService) diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp index 2960377da6..e8d5ea20e0 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp @@ -1,1161 +1,1161 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PAImageProcessing.h" // Qt #include #include #include #include //mitk image #include -#include "mitkPhotoacousticImage.h" +#include "mitkPhotoacousticFilterService.h" #include "mitkPhotoacousticBeamformingFilter.h" //other #include #include #include const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; -PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false), m_FilterBank(mitk::PhotoacousticImage::New()) +PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false), m_FilterBank(mitk::PhotoacousticFilterService::New()) { qRegisterMetaType(); qRegisterMetaType(); } void PAImageProcessing::SetFocus() { m_Controls.buttonApplyBModeFilter->setFocus(); } void PAImageProcessing::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonApplyBModeFilter, SIGNAL(clicked()), this, SLOT(StartBmodeThread())); connect(m_Controls.DoResampling, SIGNAL(clicked()), this, SLOT(UseResampling())); connect(m_Controls.Logfilter, SIGNAL(clicked()), this, SLOT(UseLogfilter())); connect(m_Controls.ResamplingValue, SIGNAL(valueChanged(double)), this, SLOT(SetResampling())); connect(m_Controls.buttonApplyBeamforming, SIGNAL(clicked()), this, SLOT(StartBeamformingThread())); connect(m_Controls.buttonApplyCropFilter, SIGNAL(clicked()), this, SLOT(StartCropThread())); connect(m_Controls.buttonApplyBandpass, SIGNAL(clicked()), this, SLOT(StartBandpassThread())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UseImageSpacing())); connect(m_Controls.ScanDepth, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBeamforming())); connect(m_Controls.BPSpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBandpass())); connect(m_Controls.Samples, SIGNAL(valueChanged(int)), this, SLOT(UpdateImageInfo())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UpdateImageInfo())); connect(m_Controls.boundLow, SIGNAL(valueChanged(int)), this, SLOT(LowerSliceBoundChanged())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpperSliceBoundChanged())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(SliceBoundsEnabled())); connect(m_Controls.BatchProcessing, SIGNAL(clicked()), this, SLOT(BatchProcessing())); connect(m_Controls.StepBeamforming, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepCropping, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBandpass, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBMode, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); UpdateSaveBoxes(); m_Controls.DoResampling->setChecked(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.progressBar->setMinimum(0); m_Controls.progressBar->setMaximum(100); m_Controls.progressBar->setVisible(false); m_Controls.UseImageSpacing->setToolTip("Image spacing of y-Axis must be in us, x-Axis in mm."); m_Controls.UseImageSpacing->setToolTipDuration(5000); m_Controls.ProgressInfo->setVisible(false); m_Controls.UseBP->hide(); m_Controls.UseGPUBmode->hide(); #ifndef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBf->setChecked(false); m_Controls.UseGPUBmode->setEnabled(false); m_Controls.UseGPUBmode->setChecked(false); #endif UseImageSpacing(); } void PAImageProcessing::ChangedSOSBandpass() { m_Controls.SpeedOfSound->setValue(m_Controls.BPSpeedOfSound->value()); } void PAImageProcessing::ChangedSOSBeamforming() { m_Controls.BPSpeedOfSound->setValue(m_Controls.SpeedOfSound->value()); } std::vector splitpath( const std::string& str , const std::set delimiters) { std::vector result; char const* pch = str.c_str(); char const* start = pch; for (; *pch; ++pch) { if (delimiters.find(*pch) != delimiters.end()) { if (start != pch) { std::string str(start, pch); result.push_back(str); } else { result.push_back(""); } start = pch + 1; } } result.push_back(start); return result; } void PAImageProcessing::UpdateSaveBoxes() { if (m_Controls.StepBeamforming->isChecked()) m_Controls.SaveBeamforming->setEnabled(true); else m_Controls.SaveBeamforming->setEnabled(false); if (m_Controls.StepCropping->isChecked()) m_Controls.SaveCropping->setEnabled(true); else m_Controls.SaveCropping->setEnabled(false); if (m_Controls.StepBandpass->isChecked()) m_Controls.SaveBandpass->setEnabled(true); else m_Controls.SaveBandpass->setEnabled(false); if (m_Controls.StepBMode->isChecked()) m_Controls.SaveBMode->setEnabled(true); else m_Controls.SaveBMode->setEnabled(false); } void PAImageProcessing::BatchProcessing() { QFileDialog LoadDialog(nullptr, "Select Files to be processed"); LoadDialog.setFileMode(QFileDialog::FileMode::ExistingFiles); LoadDialog.setNameFilter(tr("Images (*.nrrd)")); LoadDialog.setViewMode(QFileDialog::Detail); QStringList fileNames; if (LoadDialog.exec()) fileNames = LoadDialog.selectedFiles(); QString saveDir = QFileDialog::getExistingDirectory(nullptr, tr("Select Directory To Save To"), "", QFileDialog::ShowDirsOnly | QFileDialog::DontResolveSymlinks); DisableControls(); std::set delims{'/'}; bool doSteps[] = { m_Controls.StepBeamforming->isChecked(), m_Controls.StepCropping->isChecked() , m_Controls.StepBandpass->isChecked(), m_Controls.StepBMode->isChecked() }; bool saveSteps[] = { m_Controls.SaveBeamforming->isChecked(), m_Controls.SaveCropping->isChecked() , m_Controls.SaveBandpass->isChecked(), m_Controls.SaveBMode->isChecked() }; for (int fileNumber = 0; fileNumber < fileNames.size(); ++fileNumber) { m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("loading file"); QString filename = fileNames.at(fileNumber); auto split = splitpath(filename.toStdString(), delims); std::string imageName = split.at(split.size()-1); // remove ".nrrd" imageName = imageName.substr(0, imageName.size()-5); mitk::Image::Pointer image = mitk::IOUtil::Load(filename.toStdString().c_str()); UpdateBFSettings(image); // Beamforming if (doSteps[0]) { std::function progressHandle = [this](int progress, std::string progressInfo) { this->UpdateProgress(progress, progressInfo); }; m_Controls.progressBar->setValue(100); std::string errorMessage = ""; image = m_FilterBank->ApplyBeamforming(image, BFconfig, errorMessage, progressHandle); if (saveSteps[0]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " beamformed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Cropping if (doSteps[1]) { m_Controls.ProgressInfo->setText("cropping image"); image = m_FilterBank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, image->GetDimension(2) - 1); if (saveSteps[1]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " cropped" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bandpass if (doSteps[2]) { m_Controls.ProgressInfo->setText("applying bandpass"); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } image = m_FilterBank->BandpassFilter(image, recordTime, BPHighPass, BPLowPass, m_Controls.BPFalloff->value()); if (saveSteps[2]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bandpassed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bmode if (doSteps[3]) { m_Controls.ProgressInfo->setText("applying bmode filter"); bool useGPU = m_Controls.UseGPUBmode->isChecked(); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") - image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU, m_UseLogfilter, m_ResampleSpacing); + image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::Abs, useGPU, m_UseLogfilter, m_ResampleSpacing); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") - image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::EnvelopeDetection, useGPU, m_UseLogfilter, m_ResampleSpacing); + image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, useGPU, m_UseLogfilter, m_ResampleSpacing); if (saveSteps[3]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bmode" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } m_Controls.progressBar->setVisible(false); } EnableControls(); } void PAImageProcessing::StartBeamformingThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing beamforming for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("started"); m_Controls.buttonApplyBeamforming->setText("working..."); DisableControls(); BeamformingThread *thread = new BeamformingThread(); connect(thread, &BeamformingThread::result, this, &PAImageProcessing::HandleBeamformingResults); connect(thread, &BeamformingThread::updateProgress, this, &PAImageProcessing::UpdateProgress); connect(thread, &BeamformingThread::message, this, &PAImageProcessing::PAMessageBox); connect(thread, &BeamformingThread::finished, thread, &QObject::deleteLater); thread->setConfig(BFconfig); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Beamforming"; thread->start(); } } } void PAImageProcessing::HandleBeamformingResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; if (BFconfig.Algorithm == mitk::BeamformingSettings::BeamformingAlgorithm::DAS) newNodeName << "DAS bf, "; else if (BFconfig.Algorithm == mitk::BeamformingSettings::BeamformingAlgorithm::DMAS) newNodeName << "DMAS bf, "; if (BFconfig.DelayCalculationMethod == mitk::BeamformingSettings::DelayCalc::QuadApprox) newNodeName << "q. delay"; if (BFconfig.DelayCalculationMethod == mitk::BeamformingSettings::DelayCalc::Spherical) newNodeName << "s. delay"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); levelWindow.SetAuto(image, true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.ProgressInfo->setVisible(false); m_Controls.buttonApplyBeamforming->setText("Apply Beamforming"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews(image->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartBmodeThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing image processing for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyBModeFilter->setText("working..."); DisableControls(); BmodeThread *thread = new BmodeThread(); connect(thread, &BmodeThread::result, this, &PAImageProcessing::HandleBmodeResults); connect(thread, &BmodeThread::finished, thread, &QObject::deleteLater); bool useGPU = m_Controls.UseGPUBmode->isChecked(); if(m_Controls.BModeMethod->currentText() == "Absolute Filter") - thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU); + thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticFilterService::BModeMethod::Abs, useGPU); else if(m_Controls.BModeMethod->currentText() == "Envelope Detection") - thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::EnvelopeDetection, useGPU); + thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, useGPU); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Processing"; thread->start(); } } } void PAImageProcessing::HandleBmodeResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "B-Mode"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.buttonApplyBModeFilter->setText("Apply B-mode Filter"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartCropThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing image cropping for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyCropFilter->setText("working..."); DisableControls(); CropThread *thread = new CropThread(); connect(thread, &CropThread::result, this, &PAImageProcessing::HandleCropResults); connect(thread, &CropThread::finished, thread, &QObject::deleteLater); thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, image->GetDimension(2) - 1); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Cropping"; thread->start(); } } } void PAImageProcessing::HandleCropResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "Cropped"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); m_Controls.buttonApplyCropFilter->setText("Apply Crop Filter"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartBandpassThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing Bandpass filter on image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyBandpass->setText("working..."); DisableControls(); BandpassThread *thread = new BandpassThread(); connect(thread, &BandpassThread::result, this, &PAImageProcessing::HandleBandpassResults); connect(thread, &BandpassThread::finished, thread, &QObject::deleteLater); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloff->value(), recordTime); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Bandpass filter"; thread->start(); } } } void PAImageProcessing::HandleBandpassResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "Bandpassed"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); m_Controls.buttonApplyBandpass->setText("Apply Bandpass"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::SliceBoundsEnabled() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); } } void PAImageProcessing::UpperSliceBoundChanged() { if(m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundLow->setValue(m_Controls.boundHigh->value()); } } void PAImageProcessing::LowerSliceBoundChanged() { if (m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundHigh->setValue(m_Controls.boundLow->value()); } } void PAImageProcessing::UpdateProgress(int progress, std::string progressInfo) { if (progress < 100) m_Controls.progressBar->setValue(progress); else m_Controls.progressBar->setValue(100); m_Controls.ProgressInfo->setText(progressInfo.c_str()); qApp->processEvents(); } void PAImageProcessing::PAMessageBox(std::string message) { if (0 != message.compare("noMessage")) { QMessageBox msgBox; msgBox.setText(message.c_str()); msgBox.exec(); } } void PAImageProcessing::UpdateImageInfo() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { // beamforming configs if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ElementCount->setValue(image->GetDimension(0)); m_Controls.Pitch->setValue(image->GetGeometry()->GetSpacing()[0]); } m_Controls.boundLow->setMaximum(image->GetDimension(2) - 1); m_Controls.boundHigh->setMaximum(image->GetDimension(2) - 1); UpdateBFSettings(image); m_Controls.CutoffBeforeBF->setValue(0.000001 / BFconfig.TimeSpacing); // 1us standard offset for our transducer std::stringstream frequency; float maxFrequency = (1 / BFconfig.TimeSpacing) * image->GetDimension(1) / 2 / 2 / 1000; frequency << maxFrequency / 1000000; //[MHz] frequency << "MHz"; m_Controls.BPhigh->setMaximum(maxFrequency / 1000000); m_Controls.BPlow->setMaximum(maxFrequency / 1000000); frequency << " is the maximal allowed frequency for the selected image."; m_Controls.BPhigh->setToolTip(frequency.str().c_str()); m_Controls.BPlow->setToolTip(frequency.str().c_str()); m_Controls.BPhigh->setToolTipDuration(5000); m_Controls.BPlow->setToolTipDuration(5000); } } } void PAImageProcessing::OnSelectionChanged( berry::IWorkbenchPart::Pointer /*source*/, const QList& nodes ) { // iterate all selected objects, adjust warning visibility foreach( mitk::DataNode::Pointer node, nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_Controls.labelWarning->setVisible( false ); m_Controls.buttonApplyBModeFilter->setEnabled( true ); m_Controls.labelWarning2->setVisible(false); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.labelWarning3->setVisible(false); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.labelWarning4->setVisible(false); m_Controls.buttonApplyBeamforming->setEnabled(true); UpdateImageInfo(); return; } } m_Controls.labelWarning->setVisible( true ); m_Controls.buttonApplyBModeFilter->setEnabled( false ); m_Controls.labelWarning2->setVisible(true); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.labelWarning3->setVisible(true); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.labelWarning4->setVisible(true); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseResampling() { if (m_Controls.DoResampling->isChecked()) { m_Controls.ResamplingValue->setEnabled(true); m_ResampleSpacing = m_Controls.ResamplingValue->value(); } else { m_Controls.ResamplingValue->setEnabled(false); m_ResampleSpacing = 0; } } void PAImageProcessing::UseLogfilter() { m_UseLogfilter = m_Controls.Logfilter->isChecked(); } void PAImageProcessing::SetResampling() { m_ResampleSpacing = m_Controls.ResamplingValue->value(); } void PAImageProcessing::UpdateBFSettings(mitk::Image::Pointer image) { if ("DAS" == m_Controls.BFAlgorithm->currentText()) BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; else if ("sDMAS" == m_Controls.BFAlgorithm->currentText()) BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; if ("Quad. Approx." == m_Controls.DelayCalculation->currentText()) { BFconfig.DelayCalculationMethod = mitk::BeamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Wave" == m_Controls.DelayCalculation->currentText()) { BFconfig.DelayCalculationMethod = mitk::BeamformingSettings::DelayCalc::Spherical; } if ("Von Hann" == m_Controls.Apodization->currentText()) { BFconfig.Apod = mitk::BeamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { BFconfig.Apod = mitk::BeamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { BFconfig.Apod = mitk::BeamformingSettings::Apodization::Box; } BFconfig.Pitch = m_Controls.Pitch->value() / 1000; // [m] BFconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] BFconfig.SamplesPerLine = m_Controls.Samples->value(); BFconfig.ReconstructionLines = m_Controls.Lines->value(); BFconfig.TransducerElements = m_Controls.ElementCount->value(); BFconfig.apodizationArraySize = m_Controls.Lines->value(); BFconfig.Angle = m_Controls.Angle->value(); // [deg] BFconfig.UseBP = m_Controls.UseBP->isChecked(); BFconfig.UseGPU = m_Controls.UseGPUBf->isChecked(); BFconfig.upperCutoff = m_Controls.CutoffBeforeBF->value(); if (m_Controls.UseImageSpacing->isChecked()) { BFconfig.RecordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000; // [s] BFconfig.TimeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000; MITK_INFO << "Calculated Scan Depth of " << BFconfig.RecordTime * BFconfig.SpeedOfSound * 100 / 2 << "cm"; } else { BFconfig.RecordTime = 2 * m_Controls.ScanDepth->value() / 1000 / BFconfig.SpeedOfSound; // [s] BFconfig.TimeSpacing = BFconfig.RecordTime / image->GetDimension(1); } if ("US Image" == m_Controls.ImageType->currentText()) { BFconfig.isPhotoacousticImage = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { BFconfig.isPhotoacousticImage = true; } BFconfig.partial = m_Controls.Partial->isChecked(); BFconfig.CropBounds[0] = m_Controls.boundLow->value(); BFconfig.CropBounds[1] = m_Controls.boundHigh->value(); } void PAImageProcessing::EnableControls() { m_Controls.BatchProcessing->setEnabled(true); m_Controls.StepBeamforming->setEnabled(true); m_Controls.StepBandpass->setEnabled(true); m_Controls.StepCropping->setEnabled(true); m_Controls.StepBMode->setEnabled(true); UpdateSaveBoxes(); m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); m_Controls.BModeMethod->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.CutoffAbove->setEnabled(true); m_Controls.CutoffBelow->setEnabled(true); m_Controls.CutoffBeforeBF->setEnabled(true); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.BPSpeedOfSound->setEnabled(true); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.Partial->setEnabled(true); m_Controls.boundHigh->setEnabled(true); m_Controls.boundLow->setEnabled(true); m_Controls.BFAlgorithm->setEnabled(true); m_Controls.DelayCalculation->setEnabled(true); m_Controls.ImageType->setEnabled(true); m_Controls.Apodization->setEnabled(true); m_Controls.UseBP->setEnabled(true); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(true); m_Controls.UseGPUBmode->setEnabled(true); #endif m_Controls.BPhigh->setEnabled(true); m_Controls.BPlow->setEnabled(true); m_Controls.BPFalloff->setEnabled(true); m_Controls.UseImageSpacing->setEnabled(true); UseImageSpacing(); m_Controls.Pitch->setEnabled(true); m_Controls.ElementCount->setEnabled(true); m_Controls.SpeedOfSound->setEnabled(true); m_Controls.Samples->setEnabled(true); m_Controls.Lines->setEnabled(true); m_Controls.Angle->setEnabled(true); m_Controls.buttonApplyBeamforming->setEnabled(true); } void PAImageProcessing::DisableControls() { m_Controls.BatchProcessing->setEnabled(false); m_Controls.StepBeamforming->setEnabled(false); m_Controls.StepBandpass->setEnabled(false); m_Controls.StepCropping->setEnabled(false); m_Controls.StepBMode->setEnabled(false); m_Controls.SaveBeamforming->setEnabled(false); m_Controls.SaveBandpass->setEnabled(false); m_Controls.SaveCropping->setEnabled(false); m_Controls.SaveBMode->setEnabled(false); m_Controls.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); m_Controls.BModeMethod->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.CutoffAbove->setEnabled(false); m_Controls.CutoffBelow->setEnabled(false); m_Controls.CutoffBeforeBF->setEnabled(false); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.BPSpeedOfSound->setEnabled(false); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.Partial->setEnabled(false); m_Controls.boundHigh->setEnabled(false); m_Controls.boundLow->setEnabled(false); m_Controls.BFAlgorithm->setEnabled(false); m_Controls.DelayCalculation->setEnabled(false); m_Controls.ImageType->setEnabled(false); m_Controls.Apodization->setEnabled(false); m_Controls.UseBP->setEnabled(false); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBmode->setEnabled(false); #endif m_Controls.BPhigh->setEnabled(false); m_Controls.BPlow->setEnabled(false); m_Controls.BPFalloff->setEnabled(false); m_Controls.UseImageSpacing->setEnabled(false); m_Controls.ScanDepth->setEnabled(false); m_Controls.Pitch->setEnabled(false); m_Controls.ElementCount->setEnabled(false); m_Controls.SpeedOfSound->setEnabled(false); m_Controls.Samples->setEnabled(false); m_Controls.Lines->setEnabled(false); m_Controls.Angle->setEnabled(false); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseImageSpacing() { if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ScanDepth->setDisabled(true); } else { m_Controls.ScanDepth->setEnabled(true); } } #include void BeamformingThread::run() { mitk::Image::Pointer resultImage = mitk::Image::New(); mitk::Image::Pointer resultImageBuffer; std::string errorMessage = ""; std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; resultImageBuffer = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, errorMessage, progressHandle); mitk::ImageReadAccessor copy(resultImageBuffer); resultImage->Initialize(resultImageBuffer); resultImage->SetSpacing(resultImageBuffer->GetGeometry()->GetSpacing()); resultImage->SetImportVolume(const_cast(copy.GetData()), 0, 0, mitk::Image::CopyMemory); emit result(resultImage); emit message(errorMessage); } void BeamformingThread::setConfig(mitk::BeamformingSettings BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BmodeThread::run() { mitk::Image::Pointer resultImage; resultImage = m_FilterBank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseGPU, m_UseLogfilter, m_ResampleSpacing); emit result(resultImage); } -void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method, bool useGPU) +void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticFilterService::BModeMethod method, bool useGPU) { m_UseLogfilter = useLogfilter; m_ResampleSpacing = resampleSpacing; m_Method = method; m_UseGPU = useGPU; } void BmodeThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void CropThread::run() { mitk::Image::Pointer resultImage; resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, m_CutSliceFirst, m_CutSliceLast); emit result(resultImage); } void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutSliceFirst, unsigned int CutSliceLast) { m_CutAbove = CutAbove; m_CutBelow = CutBelow; m_CutSliceLast = CutSliceLast; m_CutSliceFirst = CutSliceFirst; } void CropThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BandpassThread::run() { mitk::Image::Pointer resultImage; resultImage = m_FilterBank->BandpassFilter(m_InputImage, m_RecordTime, m_BPHighPass, m_BPLowPass, m_TukeyAlpha); emit result(resultImage); } void BandpassThread::setConfig(float BPHighPass, float BPLowPass, float TukeyAlpha, float recordTime) { m_BPHighPass = BPHighPass; m_BPLowPass = BPLowPass; m_TukeyAlpha = TukeyAlpha; m_RecordTime = recordTime; } void BandpassThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h index d20b78b3ae..28cc887194 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h @@ -1,254 +1,254 @@ /*=================================================================== 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. ===================================================================*/ #ifndef PAImageProcessing_h #define PAImageProcessing_h -#include +#include #include #include #include #include "ui_PAImageProcessingControls.h" #include "mitkPhotoacousticBeamformingFilter.h" #include "mitkPhotoacousticBeamformingSettings.h" Q_DECLARE_METATYPE(mitk::Image::Pointer) Q_DECLARE_METATYPE(std::string) /*! * \brief Plugin implementing an interface for the Photoacoustic Algorithms Module * * Beamforming, Image processing as B-Mode filtering, cropping, resampling, as well as batch processing can be performed using this plugin. */ class PAImageProcessing : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; PAImageProcessing(); protected slots: void UpperSliceBoundChanged(); void LowerSliceBoundChanged(); void SliceBoundsEnabled(); void UseResampling(); void UseLogfilter(); void SetResampling(); void UseImageSpacing(); void UpdateImageInfo(); /** \brief Method called when the beamforming thread finishes; * it adds the image to a new data node and registers it to the worbench's data storage */ void HandleBeamformingResults(mitk::Image::Pointer image); /** \brief Beamforming is being performed in a separate thread to keep the workbench from freezing. */ void StartBeamformingThread(); /** \brief Method called when the B-mode filter thread finishes; * it adds the image to a new data node and registers it to the worbench's data storage */ void HandleBmodeResults(mitk::Image::Pointer image); /** \brief B-mode filtering is being performed in a separate thread to keep the workbench from freezing. */ void StartBmodeThread(); /** \brief Method called when the Cropping thread finishes; * it adds the image to a new data node and registers it to the worbench's data storage */ void HandleCropResults(mitk::Image::Pointer image); /** \brief Cropping is being performed in a separate thread to keep the workbench from freezing. */ void StartCropThread(); /** \brief Method called when the bandpass thread finishes; * it adds the image to a new data node and registers it to the worbench's data storage */ void HandleBandpassResults(mitk::Image::Pointer image); /** \brief Bandpassing is being performed in a separate thread to keep the workbench from freezing. */ void StartBandpassThread(); void UpdateProgress(int progress, std::string progressInfo); void PAMessageBox(std::string message); void BatchProcessing(); void UpdateSaveBoxes(); void ChangedSOSBandpass(); void ChangedSOSBeamforming(); protected: virtual void CreateQtPartControl(QWidget *parent) override; virtual void SetFocus() override; /** \brief called by QmitkFunctionality when DataManager's selection has changed. * On a change some parameters are internally updated to calculate bounds for GUI elements as the slice selector for beamforming or * the bandpass filter settings. */ virtual void OnSelectionChanged( berry::IWorkbenchPart::Pointer source, const QList& nodes ) override; /** \brief Instance of the GUI controls */ Ui::PAImageProcessingControls m_Controls; float m_ResampleSpacing; bool m_UseLogfilter; std::string m_OldNodeName; /** \brief The settings set which is used for beamforming, updated through this class. */ mitk::BeamformingSettings BFconfig; /** \brief Method for updating the BFconfig by using a selected image and the GUI configuration. */ void UpdateBFSettings(mitk::Image::Pointer image); void EnableControls(); void DisableControls(); /** \brief Class through which the filters are called. */ - mitk::PhotoacousticImage::Pointer m_FilterBank; + mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class BeamformingThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); void updateProgress(int, std::string); void message(std::string); public: void setConfig(mitk::BeamformingSettings BFconfig); void setInputImage(mitk::Image::Pointer image); - void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::BeamformingSettings m_BFconfig; mitk::Image::Pointer m_InputImage; int m_Cutoff; - mitk::PhotoacousticImage::Pointer m_FilterBank; + mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class BmodeThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: enum BModeMethod { ShapeDetection, Abs }; - void setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method, bool useGPU); + void setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticFilterService::BModeMethod method, bool useGPU); void setInputImage(mitk::Image::Pointer image); - void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::Image::Pointer m_InputImage; - mitk::PhotoacousticImage::BModeMethod m_Method; + mitk::PhotoacousticFilterService::BModeMethod m_Method; bool m_UseLogfilter; double m_ResampleSpacing; bool m_UseGPU; - mitk::PhotoacousticImage::Pointer m_FilterBank; + mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class CropThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutSliceFirst, unsigned int CutSliceLast); void setInputImage(mitk::Image::Pointer image); - void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::Image::Pointer m_InputImage; unsigned int m_CutAbove; unsigned int m_CutBelow; unsigned int m_CutSliceLast; unsigned int m_CutSliceFirst; - mitk::PhotoacousticImage::Pointer m_FilterBank; + mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class BandpassThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(float BPHighPass, float BPLowPass, float TukeyAlpha, float recordTime); void setInputImage(mitk::Image::Pointer image); - void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::Image::Pointer m_InputImage; float m_BPHighPass; float m_BPLowPass; float m_TukeyAlpha; float m_RecordTime; - mitk::PhotoacousticImage::Pointer m_FilterBank; + mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; #endif // PAImageProcessing_h