diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index 80c89f8144..e25901ea3d 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,30 +1,31 @@ file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*") set(CPP_FILES source/filters/mitkBeamformingFilter.cpp source/filters/mitkBeamformingSettings.cpp source/filters/mitkImageSliceSelectionFilter.cpp source/filters/mitkCastToFloatImageFilter.cpp source/filters/mitkCropImageFilter.cpp + source/filters/mitkBandpassFilter.cpp source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp source/utils/mitkPhotoacousticFilterService.cpp source/utils/mitkBeamformingUtils.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/mitkBandpassFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkBandpassFilter.h new file mode 100644 index 0000000000..6bf5357070 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBandpassFilter.h @@ -0,0 +1,57 @@ +/*=================================================================== + +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_BANDPASS_FILTER +#define MITK_BANDPASS_FILTER + +#include "mitkImageToImageFilter.h" +#include "MitkPhotoacousticsAlgorithmsExports.h" +#include "mitkPhotoacousticFilterService.h" + +namespace mitk { + /*! + * \brief Class implementing an mitk::ImageToImageFilter for casting any mitk image to a float image + */ + + class MITKPHOTOACOUSTICSALGORITHMS_EXPORT BandpassFilter : public ImageToImageFilter + { + public: + mitkClassMacro(BandpassFilter, ImageToImageFilter); + + itkFactorylessNewMacro(Self); + itkCloneMacro(Self); + + itkSetMacro(HighPass, float); + itkSetMacro(LowPass, float); + itkSetMacro(HighPassAlpha, float); + itkSetMacro(LowPassAlpha, float); + + protected: + BandpassFilter(); + + ~BandpassFilter() override; + + float m_HighPass; + float m_LowPass; + float m_HighPassAlpha; + float m_LowPassAlpha; + mitk::PhotoacousticFilterService::Pointer m_FilterService; + + void GenerateData() override; + }; +} // namespace mitk + +#endif //MITK_CAST_TO_FLOAT_IMAGE_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h index fa5d0b81be..cbe872a42d 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h @@ -1,131 +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 mitkPhotoacousticFilterService_H_HEADER_INCLUDED #define mitkPhotoacousticFilterService_H_HEADER_INCLUDED #include "itkObject.h" #include "mitkCommon.h" #include "mitkImage.h" #include #include "mitkBeamformingSettings.h" #include "mitkBeamformingFilter.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 PhotoacousticFilterService : public itk::Object { public: 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 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 UseLogFilter = false); /** \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, double outputSpacing[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 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::Pointer config, 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 alphaHighPass The high pass 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. * @param alphaLowPass The low passtukey 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, + mitk::Image::Pointer ApplyBandpassFilter(mitk::Image::Pointer data, float BPHighPass, float BPLowPass, float alphaHighPass, float alphaLowPass); protected: 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, - float cutoffFrequencyPixelHighPass, - float cutoffFrequencyPixelLowPass, - float alphaHighPass, float alphaLowPass); + mitk::Image::Pointer ConvertToFloat(mitk::Image::Pointer); }; } // namespace mitk #endif /* mitkPhotoacousticFilterService_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp new file mode 100644 index 0000000000..aaed5facfd --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp @@ -0,0 +1,237 @@ +/*=================================================================== +mitkCastToFloatImageFilter +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 "mitkBandpassFilter.h" + +#include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" +#include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" +#include "mitkImageCast.h" +#include "itkMultiplyImageFilter.h" +#include "mitkITKImageImport.h" +#include "mitkIOUtil.h" + +mitk::BandpassFilter::BandpassFilter() : + m_HighPass(0), + m_LowPass(50), + m_HighPassAlpha(1), + m_LowPassAlpha(1), + m_FilterService(mitk::PhotoacousticFilterService::New()) +{ + MITK_INFO << "Instantiating BandpassFilter..."; + SetNumberOfIndexedInputs(1); + SetNumberOfIndexedOutputs(1); + MITK_INFO << "Instantiating BandpassFilter...[Done]"; +} + +mitk::BandpassFilter::~BandpassFilter() +{ + MITK_INFO << "Destructed BandpassFilter."; +} + +itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, + float cutoffFrequencyPixelHighPass, + float cutoffFrequencyPixelLowPass, + float alphaHighPass, float alphaLowPass) +{ + float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; + + float width = cutoffFrequencyPixelLowPass - cutoffFrequencyPixelHighPass; + float center = cutoffFrequencyPixelHighPass / 2.0 + width / 2.0; + + for (unsigned int n = 0; n < reference->GetDimension(1); ++n) + { + imageData[reference->GetDimension(0)*n] = 0; + } + for (int n = 0; n < width; ++n) + { + imageData[reference->GetDimension(0)*n] = 1; + if (n <= (alphaHighPass*(width - 1)) / 2.0) + { + if (alphaHighPass > 0.00001) + { + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = + (1 + cos(itk::Math::pi*(2 * n / (alphaHighPass*(width - 1)) - 1))) / 2; + } + else + { + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; + } + } + else if (n >= (width - 1)*(1 - alphaLowPass / 2)) + { + if (alphaLowPass > 0.00001) + { + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = + (1 + cos(itk::Math::pi*(2 * n / (alphaLowPass*(width - 1)) + 1 - 2 / alphaLowPass))) / 2; + } + else + { + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; + } + } + } + + 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)]; + } + + 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; +} + +void mitk::BandpassFilter::GenerateData() +{ + auto input = GetInput(); + + if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) + { + MITK_ERROR << "Pixel type is not float, abort"; + mitkThrow() << "Pixel type is not float, abort"; + } + + auto output = GetOutput(); + mitk::Image::Pointer resampledInput; + + bool powerOfTwo = false; + int finalPower = 0; + for (int i = 1; pow(2, i) <= input->GetDimension(1); ++i) + { + finalPower = i; + if (pow(2, i) == input->GetDimension(1)) + { + powerOfTwo = true; + } + } + if (!powerOfTwo) + { + double spacing_x = input->GetGeometry()->GetSpacing()[0]; + double spacing_y = input->GetGeometry()->GetSpacing()[1] * + (((double)input->GetDimensions()[1]) / pow(2, finalPower + 1)); + double dim[2] = { spacing_x, spacing_y }; + resampledInput = m_FilterService->ApplyResampling(input, dim); + } + + // do a fourier transform, multiply with an appropriate window for the filter, and transform back + typedef itk::Image< float, 3 > RealImageType; + RealImageType::Pointer image; + + mitk::CastToItkImage(resampledInput, 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_ERROR << "Bandpass could not be applied"; + mitkThrow() << "bandpass error"; + } + + MITK_INFO << "HighPass: " << m_HighPass; + MITK_INFO << "LowPass: " << m_LowPass; + + float singleVoxel = 1.0f / ((float)resampledInput->GetGeometry()->GetSpacing()[1] / 1000); + float cutoffPixelHighPass = std::min((m_HighPass / singleVoxel) / 2, (float)resampledInput->GetDimension(1) / 2.0f); + float cutoffPixelLowPass = std::min((m_LowPass / singleVoxel) / 2, (float)resampledInput->GetDimension(1) / 2.0f - cutoffPixelHighPass); + + MITK_INFO << "cutoffPixelHighPass: " << cutoffPixelHighPass; + MITK_INFO << "cutoffPixelLowPass: " << cutoffPixelLowPass; + + RealImageType::Pointer fftMultiplicator = BPFunction(resampledInput, + cutoffPixelHighPass, + cutoffPixelLowPass, + m_HighPassAlpha, + m_LowPassAlpha); + + mitk::Image::Pointer testImage; + mitk::CastToMitkImage(fftMultiplicator, testImage); + mitk::IOUtil::Save(testImage, "G:/fftImage.nrrd"); + + typedef itk::MultiplyImageFilter< ComplexImageType, + RealImageType, + ComplexImageType > + MultiplyFilterType; + MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); + multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); + multiplyFilter->SetInput2(fftMultiplicator); + + typedef itk::FFT1DComplexConjugateToRealImageFilter< ComplexImageType, RealImageType > InverseFilterType; + InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New(); + inverseFFTFilter->SetInput(multiplyFilter->GetOutput()); + inverseFFTFilter->SetDirection(1); + + GrabItkImageMemory(inverseFFTFilter->GetOutput(), output); + mitk::IOUtil::Save(output, "G:/bandpassedImage.nrrd"); +} diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp index c53ea9da18..31b30a8f50 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp @@ -1,298 +1,298 @@ /*=================================================================== mitkBeamformingFilter 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 "mitkBeamformingFilter.h" #include "mitkBeamformingUtils.h" mitk::BeamformingFilter::BeamformingFilter(mitk::BeamformingSettings::Pointer settings) : m_OutputData(nullptr), m_InputData(nullptr), m_Conf(settings) { MITK_INFO << "Instantiating BeamformingFilter..."; this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(m_Conf); MITK_INFO << "Instantiating BeamformingFilter...[Done]"; } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { MITK_INFO << "Destructed BeamformingFilter"; } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); } 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; unsigned int dim[] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements() * 1000 / m_Conf->GetReconstructionLines(); spacing[1] = (m_Conf->GetTimeSpacing() * m_Conf->GetInputDim()[1]) / 2 * m_Conf->GetSpeedOfSound() * 1000 / m_Conf->GetSamplesPerLine(); spacing[2] = 1; output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::BeamformingFilter::GenerateData() { mitk::Image::Pointer input = this->GetInput(); if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type of input needs to be float for this filter to work."; mitkThrow() << "Pixel type of input needs to be float for this filter to work."; } GenerateOutputInformation(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf->GetUseGPU()) { 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)); m_InputData = (float*)inputReadAccessor.GetData(); m_OutputData = new float[m_Conf->GetReconstructionLines()*m_Conf->GetSamplesPerLine()]; // 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->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::sDMAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::sDMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::sDMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, 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; + mitkThrow() << "Pixel type is not float, abort"; } 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->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 4) // single output image (float) > availableMemory - (unsigned long)(m_Conf->GetReconstructionLines() / 2 * m_Conf->GetSamplesPerLine()) * 2 - // Delays buffer (unsigned short) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 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(); unsigned int num_Slices = 1; if (i == batches - 1 && (input->GetDimension(2) % batchSize > 0)) { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDimLast); num_Slices = batchDimLast[2]; } else { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDim); num_Slices = batchDim[2]; } inputBatch->SetSpacing(input->GetGeometry()->GetSpacing()); inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i])); m_BeamformingOclFilter->SetApodisation(m_Conf->GetApodizationFunction(), m_Conf->GetApodizationArraySize()); m_BeamformingOclFilter->SetInput(inputBatch); m_BeamformingOclFilter->Update(); void* out = m_BeamformingOclFilter->GetOutput(); for (unsigned int slice = 0; slice < num_Slices; ++slice) { output->SetImportSlice( &(((float*)out)[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * 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->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * m_Conf->GetInputDim()[2]]; output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } #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; } diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp index 8c2207663f..1a86a6c7a4 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp @@ -1,364 +1,210 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPhotoacousticFilterService.h" -#include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" -#include "../ITKFilter/itkPhotoacousticBModeImageFilter.h" -#include "mitkImageCast.h" + #include "mitkITKImageImport.h" -#include "mitkBeamformingFilter.h" #include -#include #include #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" #include "mitkConvert2Dto3DImageFilter.h" #include +#include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" +#include "../ITKFilter/itkPhotoacousticBModeImageFilter.h" +#include // 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" +#include "mitkImageCast.h" mitk::PhotoacousticFilterService::PhotoacousticFilterService() { MITK_INFO << "[PhotoacousticFilterService] created filter service"; } mitk::PhotoacousticFilterService::~PhotoacousticFilterService() { MITK_INFO << "[PhotoacousticFilterService] destructed filter service"; } -mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseLogFilter) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBmodeFilter( + mitk::Image::Pointer inputImage, + BModeMethod method, bool UseLogFilter) { // 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; - if (inputImage.IsNull() || - !(inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" - || inputImage->GetPixelType().GetTypeAsString() == " (float)")) - { - MITK_ERROR << "BMode Filter can only handle float image types."; - mitkThrow() << "BMode Filter can only handle float image types."; - } + auto floatImage = ConvertToFloat(inputImage); if (method == BModeMethod::Abs) { PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); filter->UseLogFilter(UseLogFilter); - filter->SetInput(inputImage); + filter->SetInput(floatImage); filter->Update(); return filter->GetOutput(); } 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); + mitk::CastToItkImage(floatImage, itkImage); if (UseLogFilter) { bModeFilter->SetInput(itkImage); bModeFilter->SetDirection(1); itkImage = bModeFilter->GetOutput(); } else { photoacousticBModeFilter->SetInput(itkImage); photoacousticBModeFilter->SetDirection(1); itkImage = photoacousticBModeFilter->GetOutput(); } return mitk::GrabItkImageMemory(itkImage); } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyResampling( mitk::Image::Pointer inputImage, double outputSpacing[2]) { typedef itk::Image< float, 3 > itkFloatImageType; + auto floatImage = ConvertToFloat(inputImage); + typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; - mitk::CastToItkImage(inputImage, itkImage); + mitk::CastToItkImage(floatImage, itkImage); itkFloatImageType::SpacingType outputSpacingItk; itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSizeItk = inputSizeItk; outputSpacingItk[0] = outputSpacing[0]; outputSpacingItk[1] = outputSpacing[1]; outputSpacingItk[2] = itkImage->GetSpacing()[2]; - outputSizeItk[0] = outputSizeItk[0] * (inputImage->GetGeometry()->GetSpacing()[0] / outputSpacing[0]); - outputSizeItk[1] = outputSizeItk[1] * (inputImage->GetGeometry()->GetSpacing()[1] / outputSpacing[1]); + outputSizeItk[0] = outputSizeItk[0] * (floatImage->GetGeometry()->GetSpacing()[0] / outputSpacing[0]); + outputSizeItk[1] = outputSizeItk[1] * (floatImage->GetGeometry()->GetSpacing()[1] / outputSpacing[1]); resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSizeItk); resampleImageFilter->SetOutputSpacing(outputSpacingItk); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyCropping( mitk::Image::Pointer inputImage, int above, int below, int right, int left, int zStart, int zEnd) { + auto floatImage = ConvertToFloat(inputImage); mitk::CropImageFilter::Pointer cropImageFilter = mitk::CropImageFilter::New(); - cropImageFilter->SetInput(inputImage); + cropImageFilter->SetInput(floatImage); cropImageFilter->SetXPixelsCropStart(left); cropImageFilter->SetXPixelsCropEnd(right); cropImageFilter->SetYPixelsCropStart(above); cropImageFilter->SetYPixelsCropEnd(below); cropImageFilter->SetZPixelsCropStart(zStart); cropImageFilter->SetZPixelsCropEnd(zEnd); cropImageFilter->Update(); return cropImageFilter->GetOutput(); } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBeamforming( mitk::Image::Pointer inputImage, BeamformingSettings::Pointer config, std::function progressHandle) { Image::Pointer processedImage = mitk::Image::New(); if (inputImage->GetDimension() != 3) { mitk::Convert2Dto3DImageFilter::Pointer dimensionImageFilter = mitk::Convert2Dto3DImageFilter::New(); dimensionImageFilter->SetInput(inputImage); dimensionImageFilter->Update(); processedImage = dimensionImageFilter->GetOutput(); } else { processedImage = inputImage; } - // perform the beamforming m_BeamformingFilter = mitk::BeamformingFilter::New(config); - m_BeamformingFilter->SetInput(processedImage); + m_BeamformingFilter->SetInput(ConvertToFloat(processedImage)); m_BeamformingFilter->SetProgressHandle(progressHandle); m_BeamformingFilter->UpdateLargestPossibleRegion(); processedImage = m_BeamformingFilter->GetOutput(); return processedImage; } -mitk::Image::Pointer mitk::PhotoacousticFilterService::BandpassFilter( - mitk::Image::Pointer data, float recordTime, +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBandpassFilter( + mitk::Image::Pointer data, float BPHighPass, float BPLowPass, float alphaHighPass, float alphaLowPass) { - bool powerOfTwo = false; - int finalPower = 0; - for (int i = 1; pow(2, i) <= data->GetDimension(1); ++i) - { - finalPower = i; - if (pow(2, i) == data->GetDimension(1)) - { - powerOfTwo = true; - } - } - if (!powerOfTwo) - { - double spacing_x = data->GetGeometry()->GetSpacing()[0]; - double spacing_y = data->GetGeometry()->GetSpacing()[1] * (((double)data->GetDimensions()[1]) / pow(2, finalPower + 1)); - double dim[2] = { spacing_x, spacing_y }; - data = ApplyResampling(data, dim); - } - - MITK_INFO << data->GetDimension(0); - MITK_INFO << data->GetDimension(1); - - // do a fourier transform, multiply with an appropriate window for the filter, and transform back - typedef itk::Image< float, 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_ERROR << "Bandpass could not be applied"; - mitkThrow() << "bandpass error"; - } - - float singleVoxel = 1 / (recordTime / data->GetDimension(1)) / 2 / 1000; - float cutoffPixelHighPass = std::min(BPHighPass / singleVoxel, (float)data->GetDimension(1) / 2); - float cutoffPixelLowPass = std::min(BPLowPass / singleVoxel, (float)data->GetDimension(1) / 2 - cutoffPixelHighPass); - - RealImageType::Pointer fftMultiplicator = BPFunction(data, cutoffPixelHighPass, cutoffPixelLowPass, alphaHighPass, alphaLowPass); - - typedef itk::MultiplyImageFilter< ComplexImageType, - RealImageType, - ComplexImageType > - MultiplyFilterType; - MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); - multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); - multiplyFilter->SetInput2(fftMultiplicator); - - typedef itk::FFT1DComplexConjugateToRealImageFilter< ComplexImageType, RealImageType > InverseFilterType; - InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New(); - inverseFFTFilter->SetInput(multiplyFilter->GetOutput()); - inverseFFTFilter->SetDirection(1); - - return GrabItkImageMemory(inverseFFTFilter->GetOutput()); + auto floatData = ConvertToFloat(data); + mitk::BandpassFilter::Pointer bandpassFilter = mitk::BandpassFilter::New(); + bandpassFilter->SetInput(floatData); + bandpassFilter->SetHighPass(BPHighPass); + bandpassFilter->SetLowPass(BPLowPass); + bandpassFilter->SetHighPassAlpha(alphaHighPass); + bandpassFilter->SetLowPassAlpha(alphaLowPass); + bandpassFilter->Update(); + return bandpassFilter->GetOutput(); } -itk::Image::Pointer mitk::PhotoacousticFilterService::BPFunction(mitk::Image::Pointer reference, - float cutoffFrequencyPixelHighPass, - float cutoffFrequencyPixelLowPass, - float alphaHighPass, float alphaLowPass) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ConvertToFloat(mitk::Image::Pointer inputImage) { - float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; - - float width = reference->GetDimension(1) / 2.0 - cutoffFrequencyPixelHighPass - cutoffFrequencyPixelLowPass; - float center = cutoffFrequencyPixelHighPass / 2.0 + width / 2.0; - - for (unsigned int n = 0; n < reference->GetDimension(1); ++n) - { - imageData[reference->GetDimension(0)*n] = 0; - } - for (int n = 0; n < width; ++n) - { - imageData[reference->GetDimension(0)*n] = 1; - if (n <= (alphaHighPass*(width - 1)) / 2.0) - { - if (alphaHighPass > 0.00001) - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = - (1 + cos(itk::Math::pi*(2 * n / (alphaHighPass*(width - 1)) - 1))) / 2; - } - else - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; - } - } - else if (n >= (width - 1)*(1 - alphaLowPass / 2)) - { - if (alphaLowPass > 0.00001) - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = - (1 + cos(itk::Math::pi*(2 * n / (alphaLowPass*(width - 1)) + 1 - 2 / alphaLowPass))) / 2; - } - else - { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; - } - } - } - - 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)]; - } - - for (unsigned int line = 1; line < reference->GetDimension(0); ++line) + if ((inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || + inputImage->GetPixelType().GetTypeAsString() == " (float)")) { - for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample) - { - imageData[reference->GetDimension(0)*sample + line] = imageData[reference->GetDimension(0)*sample]; - } + return inputImage; } - 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; + mitk::CastToFloatImageFilter::Pointer castToFloatImageFilter = mitk::CastToFloatImageFilter::New(); + castToFloatImageFilter->SetInput(inputImage); + castToFloatImageFilter->Update(); + return castToFloatImageFilter->GetOutput(); } 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 41a0274cac..1164580488 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,1065 +1,1020 @@ /*=================================================================== 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 "mitkPhotoacousticFilterService.h" #include "mitkCastToFloatImageFilter.h" #include "mitkBeamformingFilter.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::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()); auto BFconfig = CreateBeamformingSettings(image); // Beamforming if (doSteps[0]) { std::function progressHandle = [this](int progress, std::string progressInfo) { this->UpdateProgress(progress, progressInfo); }; m_Controls.progressBar->setValue(100); image = m_FilterBank->ApplyBeamforming(image, BFconfig, 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, 0); 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 - 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, + image = m_FilterBank->ApplyBandpassFilter(image, BPHighPass, BPLowPass, m_Controls.BPFalloffHigh->value(), m_Controls.BPFalloffLow->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"); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::Abs, m_UseLogfilter); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, m_UseLogfilter); double desiredSpacing[2]{ image->GetGeometry()->GetSpacing()[0], m_ResampleSpacing }; image = m_FilterBank->ApplyResampling(image, desiredSpacing); 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) { auto BFconfig = CreateBeamformingSettings(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"); DisableControls(); BeamformingThread *thread = new BeamformingThread(); connect(thread, &BeamformingThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BeamformingThread::updateProgress, this, &PAImageProcessing::UpdateProgress); 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::HandleResults(mitk::Image::Pointer image, std::string nameExtension) { MITK_INFO << "Handling results..."; auto newNode = mitk::DataNode::New(); newNode->SetData(image); newNode->SetName(m_OldNodeName + nameExtension); // 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); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews(image->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); MITK_INFO << "Handling results...[Done]"; } 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) { 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(); DisableControls(); BmodeThread *thread = new BmodeThread(); connect(thread, &BmodeThread::result, this, &PAImageProcessing::HandleResults); 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::PhotoacousticFilterService::BModeMethod::Abs, useGPU); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") 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::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) { 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(); DisableControls(); CropThread *thread = new CropThread(); connect(thread, &CropThread::result, this, &PAImageProcessing::HandleResults); 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::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) { 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(); DisableControls(); BandpassThread *thread = new BandpassThread(); connect(thread, &BandpassThread::result, this, &PAImageProcessing::HandleResults); 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(); + float BPHighPass = 1000000.0f * m_Controls.BPhigh->value(); // [Now in Hz] + float BPLowPass = 1000000.0f * m_Controls.BPlow->value(); // [Now in Hz] - BPHighPass = 0; - } - if (BPHighPass > maxFrequency - BPLowPass) - { - QMessageBox Msgbox; - Msgbox.setText("HighPass higher than LowPass, disabled both."); - Msgbox.exec(); - - BPHighPass = 0; - BPLowPass = 0; - } - - thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloffLow->value(), m_Controls.BPFalloffHigh->value(), recordTime); + thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloffLow->value(), m_Controls.BPFalloffHigh->value(), 0); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Bandpass filter"; thread->start(); } } } 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); float speedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] std::stringstream frequency; float timeSpacing; if (m_Controls.UseImageSpacing->isChecked()) { timeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000.0f; MITK_INFO << "Calculated Scan Depth of " << (image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000) * speedOfSound * 100 / 2 << "cm"; } else { timeSpacing = (2 * m_Controls.ScanDepth->value() / 1000 / speedOfSound) / image->GetDimension(1); } float maxFrequency = (1 / 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(); } mitk::BeamformingSettings::Pointer PAImageProcessing::CreateBeamformingSettings(mitk::Image::Pointer image) { mitk::BeamformingSettings::BeamformingAlgorithm algorithm; if ("DAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; else if ("sDMAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; mitk::BeamformingSettings::DelayCalc delay; if ("Quad. Approx." == m_Controls.DelayCalculation->currentText()) { delay = mitk::BeamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Wave" == m_Controls.DelayCalculation->currentText()) { delay = mitk::BeamformingSettings::DelayCalc::Spherical; } mitk::BeamformingSettings::Apodization apod; if ("Von Hann" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Box; } float pitchInMeters = m_Controls.Pitch->value() / 1000; // [m] float speedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] unsigned int samplesPerLine = m_Controls.Samples->value(); unsigned int reconstructionLines = m_Controls.Lines->value(); unsigned int apodizatonArraySize = m_Controls.Lines->value(); float angle = m_Controls.Angle->value(); // [deg] bool useBandpass = m_Controls.UseBP->isChecked(); bool useGPU = m_Controls.UseGPUBf->isChecked(); float timeSpacing; if (m_Controls.UseImageSpacing->isChecked()) { timeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000.0f; MITK_INFO << "Calculated Scan Depth of " << (image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000) * speedOfSound * 100 / 2 << "cm"; } else { timeSpacing = (2 * m_Controls.ScanDepth->value() / 1000 / speedOfSound) / image->GetDimension(1); } bool isPAImage; if ("US Image" == m_Controls.ImageType->currentText()) { isPAImage = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { isPAImage = true; } bool partial = m_Controls.Partial->isChecked(); unsigned int cropBoundLow = m_Controls.boundLow->value(); unsigned int cropBoundHigh = m_Controls.boundHigh->value(); unsigned int* cropBounds = new unsigned int[2]{ cropBoundLow, cropBoundHigh }; return mitk::BeamformingSettings::New(pitchInMeters, speedOfSound, timeSpacing, angle, isPAImage, samplesPerLine, reconstructionLines, 0, partial, cropBounds, image->GetDimensions(), useGPU, delay, apod, apodizatonArraySize, algorithm, useBandpass, 50, 50); } 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.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.BPFalloffLow->setEnabled(true); m_Controls.BPFalloffHigh->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.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.BPFalloffLow->setEnabled(false); m_Controls.BPFalloffHigh->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::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; - mitk::CastToFloatImageFilter::Pointer castToFloatImageFilter = mitk::CastToFloatImageFilter::New(); - castToFloatImageFilter->SetInput(m_InputImage); - castToFloatImageFilter->Update(); - m_InputImage = castToFloatImageFilter->GetOutput(); - mitk::IOUtil::Save(m_InputImage, "G:/castedImage.nrrd"); resultImageBuffer = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, progressHandle); emit result(resultImageBuffer, "_bf"); } void BeamformingThread::setConfig(mitk::BeamformingSettings::Pointer BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BmodeThread::run() { mitk::Image::Pointer resultImage = m_FilterBank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseLogfilter); double desiredSpacing[2]{ m_InputImage->GetGeometry()->GetSpacing()[0], m_ResampleSpacing }; resultImage = m_FilterBank->ApplyResampling(resultImage, desiredSpacing); emit result(resultImage, "_bmode"); } 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, 0, 0); emit result(resultImage, "_cropped"); } 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_TukeyAlphaHighPass, m_TukeyAlphaLowPass); + mitk::Image::Pointer resultImage = m_FilterBank->ApplyBandpassFilter(m_InputImage, m_BPHighPass, m_BPLowPass, m_TukeyAlphaHighPass, m_TukeyAlphaLowPass); emit result(resultImage, "_bandpassed"); } void BandpassThread::setConfig(float BPHighPass, float BPLowPass, float TukeyAlphaHighPass, float TukeyAlphaLowPass, float recordTime) { m_BPHighPass = BPHighPass; m_BPLowPass = BPLowPass; m_TukeyAlphaHighPass = TukeyAlphaHighPass; m_TukeyAlphaLowPass = TukeyAlphaLowPass; m_RecordTime = recordTime; } void BandpassThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; }