diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/DMASspherical.cl b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/DMASspherical.cl new file mode 100644 index 0000000000..d6dde1a66e --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/DMASspherical.cl @@ -0,0 +1,40 @@ +/*=================================================================== + +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. + +===================================================================*/ + +__kernel void ckBinaryThreshold( + __read_only image3d_t dSource, // input image + __global double* dDest, // output buffer + unsigned short inputDim[2], unsigned short outputDim[2], unsigned short apodArraySize, double* apodArray // parameters +) +{ + // get thread identifier + unsigned int globalPosX = get_global_id(0); + unsigned int globalPosY = get_global_id(1); + unsigned int globalPosZ = get_global_id(2); + + // create an image sampler + const sampler_t defaultSampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST ; + + // terminate non-valid threads + if ( globalPosX < uiWidth && globalPosY < uiHeight && globalPosZ < uiDepth ) + { + result = globalPosX; + + // store the result + dDest[ globalPosZ * uiWidth * uiHeight + globalPosY * uiWidth + globalPosX ] = result; + } + +} diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp new file mode 100644 index 0000000000..47019949b1 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp @@ -0,0 +1,110 @@ +/*=================================================================== + +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 "mitkPhotoacousticOCLBeamformer.h" +#include "usServiceReference.h" +/* +mitk::PhotoacousticOCLBeamformer::PhotoacousticOCLBeamformer() +: m_PixelCalculation( NULL ) +{ + this->AddSourceFile("BinaryThresholdFilter.cl"); + this->m_FilterID = "PixelCalculation"; +} + +mitk::PhotoacousticOCLBeamformer::~PhotoacousticOCLBeamformer() +{ + if ( this->m_PixelCalculation ) + { + clReleaseKernel( m_PixelCalculation ); + } +} + +void mitk::PhotoacousticOCLBeamformer::Update() +{ + //Check if context & program available + if (!this->Initialize()) + { + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + + // clean-up also the resources + resources->InvalidateStorage(); + mitkThrow() <<"Filter is not initialized. Cannot update."; + } + else{ + // Execute + this->Execute(); + } +} + +void mitk::PhotoacousticOCLBeamformer::Execute() +{ + cl_int clErr = 0; + + try + { + this->InitExec( this->m_PixelCalculation ); + } + catch( const mitk::Exception& e) + { + MITK_ERROR << "Catched exception while initializing filter: " << e.what(); + return; + } + + // set kernel arguments + clErr = clSetKernelArg( this->m_PixelCalculation, 2, sizeof(cl_int), &(this->m_LowerThreshold) ); + clErr |= clSetKernelArg( this->m_PixelCalculation, 3, sizeof(cl_int), &(this->m_UpperThreshold) ); + clErr |= clSetKernelArg( this->m_PixelCalculation, 4, sizeof(cl_int), &(this->m_OutsideValue) ); + clErr |= clSetKernelArg( this->m_PixelCalculation, 5, sizeof(cl_int), &(this->m_InsideValue) ); + CHECK_OCL_ERR( clErr ); + + // execute the filter on a 3D NDRange + this->SetWorkingSize(1, m_outputDim[0], 1, outputDim[1], 1, outputDim[2]); + this->ExecuteKernel( m_PixelCalculation, 3); + + // signalize the GPU-side data changed + m_Output->Modified( GPU_DATA ); +} + +us::Module *mitk::PhotoacousticOCLBeamformer::GetModule() +{ + return us::GetModuleContext()->GetModule(); +} + +bool mitk::PhotoacousticOCLBeamformer::Initialize() +{ + bool buildErr = true; + cl_int clErr = 0; + + if ( OclFilter::Initialize() ) + { + this->m_PixelCalculation = clCreateKernel( this->m_ClProgram, "ckPixelCalculation", &clErr); + buildErr |= CHECK_OCL_ERR( clErr ); + } + + return (OclFilter::IsInitialized() && buildErr ); +} + +void mitk::PhotoacousticOCLBeamformer::SetInput(mitk::Image::Pointer image) +{ + if(image->GetDimension() != 3) + { + mitkThrowException(mitk::Exception) << "Input for " << this->GetNameOfClass() << + " is not 3D. The filter only supports 3D. Please change your input."; + } + OclImageToImageFilter::SetInput(image); +} +*/ \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h new file mode 100644 index 0000000000..8ae501e586 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h @@ -0,0 +1,94 @@ +/*=================================================================== + +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 _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ +#define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ + +#include "..\\..\\..\\OpenCL\\mitkOclImageToImageFilter.h" +#include +#ifdef HAHA +namespace mitk +{ +class OclImageToImageFilter; + +/** Documentation + * + * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given + threshold values. + + * + * The filter requires two threshold values ( the upper and the lower threshold ) and two image values ( inside and outside ). The resulting voxel of the segmentation image is assigned the inside value 1 if the image value is between the given thresholds and the outside value otherwise. + */ + + +class PhotoacousticOCLBeamformer : public OclImageToImageFilter, public itk::Object +{ + +public: + mitkClassMacroItkParent(PhotoacousticOCLBeamformer, itk::Object); + itkNewMacro(Self); + + /** + * @brief SetInput Set the input image. Only 3D images are supported for now. + * @param image a 3D image. + * @throw mitk::Exception if the dimesion is not 3. + */ + void SetInput(Image::Pointer image); + + /** Update the filter */ + void Update(); + + void SetOutputDim( unsigned int outputDim[3]) + { + m_OutputDim[0] = outputDim[0]; + m_OutputDim[1] = outputDim[1]; + m_OutputDim[2] = outputDim[2]; + } + +protected: + + /** Constructor */ + PhotoacousticOCLBeamformer(); + + /** Destructor */ + virtual ~PhotoacousticOCLBeamformer(); + + /** Initialize the filter */ + bool Initialize(); + + void Execute(); + + mitk::PixelType GetOutputType() + { + return mitk::MakeScalarPixelType(); + } + + int GetBytesPerElem() + { + return sizeof(double); + } + + virtual us::Module* GetModule(); + +private: + /** The OpenCL kernel for the filter */ + cl_kernel m_PixelCalculation; + + unsigned int m_outputDim[3]; +}; +} +#endif +#endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp deleted file mode 100644 index a07076df4b..0000000000 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp +++ /dev/null @@ -1,563 +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. - -===================================================================*/ - -#define _USE_MATH_DEFINES - -#include "mitkPhotoacousticBeamformingDASFilter.h" -#include "mitkProperties.h" -#include "mitkImageReadAccessor.h" -#include -#include -#include -#include -#include -#include -#include "itkFFT1DComplexConjugateToRealImageFilter.h" -#include "itkFFT1DRealToComplexConjugateImageFilter.h" -#include "mitkImageCast.h" - -// needed itk image filters -#include "mitkITKImageImport.h" -#include "itkFFTShiftImageFilter.h" -#include "itkMultiplyImageFilter.h" -#include "itkComplexToModulusImageFilter.h" -#include - - -mitk::BeamformingDASFilter::BeamformingDASFilter() : m_OutputData(nullptr), m_InputData(nullptr) -{ - this->SetNumberOfIndexedInputs(1); - this->SetNumberOfRequiredInputs(1); - - m_ProgressHandle = [](int, std::string) {}; -} - -void mitk::BeamformingDASFilter::SetProgressHandle(std::function progressHandle) -{ - m_ProgressHandle = progressHandle; -} - -mitk::BeamformingDASFilter::~BeamformingDASFilter() -{ -} - -void mitk::BeamformingDASFilter::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::BeamformingDASFilter::GenerateOutputInformation() -{ - mitk::Image::ConstPointer input = this->GetInput(); - mitk::Image::Pointer output = this->GetOutput(); - - if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) - return; - - itkDebugMacro(<< "GenerateOutputInformation()"); - - if (!m_Conf.partial) - { - m_Conf.bounds[0] = 0; - m_Conf.bounds[1] = input->GetDimension(2) - 1; - } - - unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) - m_Conf.bounds[0] - ((input->GetDimension(2) - 1) - m_Conf.bounds[1]) }; - 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.RecordTime * 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::BeamformingDASFilter::GenerateData() -{ - GenerateOutputInformation(); - mitk::Image::ConstPointer input = this->GetInput(); - mitk::Image::Pointer output = this->GetOutput(); - - if (!output->IsInitialized()) - return; - - double inputDim[2] = { input->GetDimension(0), input->GetDimension(1) / ((int)m_Conf.Photoacoustic + 1) }; - // if the photoacoustic option is used, we halve the image, as only the upper part of it contains any information - double outputDim[2] = { output->GetDimension(0), output->GetDimension(1) }; - - const int apodArraySize = m_Conf.TransducerElements * 4; // set the resolution of the apodization array - double* ApodWindow; - // calculate the appropiate apodization window - if (m_Conf.Apod == beamformingSettings::Apodization::Hann) - { - ApodWindow = VonHannFunction(apodArraySize); - } - else if (m_Conf.Apod == beamformingSettings::Apodization::Hamm) - { - ApodWindow = HammFunction(apodArraySize); - } - else - { - ApodWindow = BoxFunction(apodArraySize); - } - - int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; - // the interval at which we update the gui progress bar - - auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... - - for (int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied - { - mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i + m_Conf.bounds[0])); - - m_OutputData = new double[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; - m_InputDataPuffer = new double[input->GetDimension(0)*input->GetDimension(1)]; - - // first, we convert any data to double, which we use by default - if (input->GetPixelType().GetTypeAsString() == "scalar (double)" || input->GetPixelType().GetTypeAsString() == " (double)") - { - m_InputData = (double*)inputReadAccessor.GetData(); - } - else if (input->GetPixelType().GetTypeAsString() == "scalar (short)") - { - short* InputPuffer = (short*)inputReadAccessor.GetData(); - for (int l = 0; l < inputDim[0]; ++l) - { - for (int s = 0; s < inputDim[1]; ++s) - { - m_InputDataPuffer[l*(short)inputDim[1] + s] = (double)InputPuffer[l*(short)inputDim[1] + s]; - } - } - m_InputData = m_InputDataPuffer; - } - else if (input->GetPixelType().GetTypeAsString() == "scalar (float)") - { - float* InputPuffer = (float*)inputReadAccessor.GetData(); - for (int l = 0; l < inputDim[0]; ++l) - { - for (int s = 0; s < inputDim[1]; ++s) - { - m_InputDataPuffer[l*(short)inputDim[1] + s] = (double)InputPuffer[l*(short)inputDim[1] + s]; - } - } - m_InputData = m_InputDataPuffer; - } - else - { - MITK_INFO << "Could not determine pixel type"; - return; - } - - // 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 - for (short line = 0; line < outputDim[0]; ++line) - { - if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) - { - threads[line] = std::thread(&BeamformingDASFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); - } - else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) - { - threads[line] = std::thread(&BeamformingDASFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); - } - } - // 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) / (double)output->GetDimension(2) * 100), "performing reconstruction"); - - delete[] m_OutputData; - delete[] m_InputDataPuffer; - m_OutputData = nullptr; - m_InputData = nullptr; - } - - // apply a bandpass filter, if requested - if (m_Conf.UseBP) - { - m_ProgressHandle(100, "applying bandpass"); - mitk::Image::Pointer BP = BandpassFilter(output); - - for (int i = 0; i < output->GetDimension(2); ++i) - { - mitk::ImageReadAccessor copy(BP, BP->GetSliceData(i)); - output->SetSlice(copy.GetData(), i); - } - } - - m_TimeOfHeaderInitialization.Modified(); - - auto end = std::chrono::high_resolution_clock::now(); - MITK_INFO << "DAS Beamforming of " << output->GetDimension(2) << " Images completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; -} - -void mitk::BeamformingDASFilter::Configure(beamformingSettings settings) -{ - m_Conf = settings; -} - -double* mitk::BeamformingDASFilter::VonHannFunction(int samples) -{ - double* ApodWindow = new double[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; - } - - return ApodWindow; -} - -double* mitk::BeamformingDASFilter::HammFunction(int samples) -{ - double* ApodWindow = new double[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = 0.54 - 0.46*cos(2 * M_PI*n / (samples - 1)); - } - - return ApodWindow; -} - -double* mitk::BeamformingDASFilter::BoxFunction(int samples) -{ - double* ApodWindow = new double[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = 1; - } - - return ApodWindow; -} - -void mitk::BeamformingDASFilter::DASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) -{ - double& inputS = inputDim[1]; - double& inputL = inputDim[0]; - - double& outputS = outputDim[1]; - double& outputL = outputDim[0]; - - short AddSample = 0; - short maxLine = 0; - short minLine = 0; - double delayMultiplicator = 0; - double l_i = 0; - double s_i = 0; - - double l = 0; - double x = 0; - double root = 0; - - double part = 0.07 * inputL; - double tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); - double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; - double apod_mult = 1; - - double 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; - - 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.0); - usedLines = (maxLine - minLine); - - apod_mult = apodArraySize / (maxLine - minLine); - - delayMultiplicator = pow((inputS / (m_Conf.RecordTime*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; - 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::BeamformingDASFilter::DASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) -{ - double& inputS = inputDim[1]; - double& inputL = inputDim[0]; - - double& outputS = outputDim[1]; - double& outputL = outputDim[0]; - - short AddSample = 0; - short maxLine = 0; - short minLine = 0; - double delayMultiplicator = 0; - double l_i = 0; - double s_i = 0; - - double l = 0; - double x = 0; - double root = 0; - - double part = 0.07 * inputL; - double tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); - double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; - double apod_mult = 1; - - double mult = 0; - short usedLines = (maxLine - minLine); - - //exact delay - - l_i = (double)line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (double)sample / outputS * inputS; - - 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.0); - usedLines = (maxLine - minLine); - - apod_mult = apodArraySize / (maxLine - minLine); - - for (short l_s = minLine; l_s < maxLine; ++l_s) - { - AddSample = (int)sqrt( - pow(s_i, 2) - + - pow((inputS / (m_Conf.RecordTime*m_Conf.SpeedOfSound) * ((l_s - l_i)*m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) - ); - 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; - } -} - -mitk::Image::Pointer mitk::BeamformingDASFilter::BandpassFilter(mitk::Image::Pointer data) -{ - // do a fourier transform, multiply with an appropriate window for the filter, and transform back - typedef double 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 can not be applied after beamforming"; - return data; - } - - double singleVoxel = 1 / (m_Conf.RecordTime / data->GetDimension(1)) / 2 / 1000; - double BoundHighPass = std::min(m_Conf.BPHighPass / singleVoxel, (double)data->GetDimension(1) / 2); - double BoundLowPass = std::min(m_Conf.BPLowPass / singleVoxel, (double)data->GetDimension(1) / 2 - BoundHighPass); - - int center1 = ((-BoundLowPass - BoundHighPass + data->GetDimension(1) / 2) / 2) + BoundLowPass; - int center2 = ((-BoundLowPass - BoundHighPass + data->GetDimension(1) / 2) / 2) + BoundHighPass + data->GetDimension(1) / 2; - - int width1 = -BoundLowPass - BoundHighPass + data->GetDimension(1) / 2; - int width2 = -BoundLowPass - BoundHighPass + data->GetDimension(1) / 2; - - /*MITK_INFO << "BHP " << BoundHighPass << " BLP " << BoundLowPass << "BPLP" << m_Conf.BPLowPass; - MITK_INFO << "center1 " << center1 << " width1 " << width1; - MITK_INFO << "center2 " << center2 << " width2 " << width2;*/ //DEBUG - - RealImageType::Pointer fftMultiplicator1 = BPFunction(data, width1, center1); - RealImageType::Pointer fftMultiplicator2 = BPFunction(data, width2, center2); - - typedef itk::AddImageFilter AddImageFilterType; - AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New(); - addImageFilter->SetInput1(fftMultiplicator1); - addImageFilter->SetInput2(fftMultiplicator2); - - typedef itk::FFTShiftImageFilter< RealImageType, RealImageType > FFTShiftFilterType; - FFTShiftFilterType::Pointer fftShiftFilter = FFTShiftFilterType::New(); - fftShiftFilter->SetInput(addImageFilter->GetOutput()); - - typedef itk::MultiplyImageFilter< ComplexImageType, - RealImageType, - ComplexImageType > - MultiplyFilterType; - MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); - multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); - multiplyFilter->SetInput2(fftShiftFilter->GetOutput()); - - /*itk::ComplexToModulusImageFilter::Pointer toReal = itk::ComplexToModulusImageFilter::New(); - toReal->SetInput(multiplyFilter->GetOutput()); - return GrabItkImageMemory(addImageFilter->GetOutput());*/ //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::BeamformingDASFilter::BPFunction(mitk::Image::Pointer reference, int width, int center) -{ - double* imageData = new double[reference->GetDimension(0)*reference->GetDimension(1)]; - - for (int sample = 0; sample < reference->GetDimension(1); ++sample) - { - imageData[reference->GetDimension(0)*sample] = 0; - } - - /* // tukey window - double alpha = m_Conf.BPFalloff; - - for (int n = 0; n < width; ++n) - { - if (n <= (alpha*(width - 1)) / 2) - { - imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) - 1))) / 2; - } - else if (n >= (width - 1)*(1 - alpha / 2) && n <= (width - 1)) - { - imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; - } - else - { - imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = 1; - } - }*/ - - // Butterworth-Filter - double d = center - width / 2; - double l = center + width / 2; - - for (int n = 0; n < reference->GetDimension(1); ++n) - { - imageData[reference->GetDimension(0)*n] = 1 / (1 + pow(((double)center - (double)n) / ((double)width / 2), 2 * m_Conf.ButterworthOrder)); - } - - // copy and paste to all lines - for (int line = 1; line < reference->GetDimension(0); ++line) - { - for (int sample = 0; sample < reference->GetDimension(1); ++sample) - { - imageData[reference->GetDimension(0)*sample + line] = imageData[reference->GetDimension(0)*sample]; - } - } - - typedef itk::Image< double, 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; -} diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.h b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.h deleted file mode 100644 index 3eff031875..0000000000 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.h +++ /dev/null @@ -1,103 +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 MITK_PHOTOACOUSTICS_BEAMFORMING_DMAS_FILTER -#define MITK_PHOTOACOUSTICS_BEAMFORMING_DMAS_FILTER - -#include "mitkImageToImageFilter.h" -#include - -namespace mitk { - - //##Documentation - //## @brief - //## @ingroup Process - class BeamformingDMASFilter : public ImageToImageFilter - { - public: - mitkClassMacro(BeamformingDMASFilter, ImageToImageFilter); - - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - struct beamformingSettings - { - double Pitch = 0.0003; // [m] - double SpeedOfSound = 1540; // [m/s] - unsigned int SamplesPerLine = 2048; - unsigned int ReconstructionLines = 128; - double RecordTime = 0.00006; // [s] - unsigned int TransducerElements = 128; - bool partial = false; - unsigned int bounds[2] = { 0,0 }; - - enum DelayCalc { QuadApprox, Spherical }; - DelayCalc DelayCalculationMethod = QuadApprox; - - enum Apodization { Hamm, Hann, Box }; - Apodization Apod = Hann; - - double Angle = 10.0; - bool Photoacoustic = true; - double BPHighPass = 50; - double BPLowPass = 50; - bool UseBP = true; - unsigned int ButterworthOrder = 8; - }; - - void Configure(beamformingSettings settings); - - void SetProgressHandle(std::function progressHandle); - - protected: - - BeamformingDMASFilter(); - - ~BeamformingDMASFilter(); - - virtual void GenerateInputRequestedRegion() override; - - virtual void GenerateOutputInformation() override; - - virtual void GenerateData() override; - - //##Description - //## @brief Time when Header was last initialized - itk::TimeStamp m_TimeOfHeaderInitialization; - - std::function m_ProgressHandle; - - double* VonHannFunction(int samples); - double* HammFunction(int samples); - double* BoxFunction(int samples); - - void DMASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize); - void DMASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize); - - mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data); - itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int width, int center); - - double* m_OutputData; - double* m_InputData; - double* m_InputDataPuffer; - - beamformingSettings m_Conf; - }; - -} // namespace mitk - -#endif MITK_PHOTOACOUSTICS_BEAMFORMING_DMAS_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp similarity index 70% rename from Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp rename to Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp index 98b6e3fc44..ff1da3e379 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp @@ -1,646 +1,746 @@ /*=================================================================== - +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. ===================================================================*/ + #define _USE_MATH_DEFINES -#include "mitkPhotoacousticBeamformingDMASFilter.h" +#include "mitkPhotoacousticBeamformingFilter.h" #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include +#include +#include +#include #include "itkFFT1DComplexConjugateToRealImageFilter.h" #include "itkFFT1DRealToComplexConjugateImageFilter.h" #include "mitkImageCast.h" -#include -#include // needed itk image filters #include "mitkITKImageImport.h" #include "itkFFTShiftImageFilter.h" #include "itkMultiplyImageFilter.h" #include "itkComplexToModulusImageFilter.h" #include -mitk::BeamformingDMASFilter::BeamformingDMASFilter() : m_OutputData(nullptr), m_InputData(nullptr) + +mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; } -void mitk::BeamformingDMASFilter::SetProgressHandle(std::function progressHandle) +void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } -mitk::BeamformingDMASFilter::~BeamformingDMASFilter() +mitk::BeamformingFilter::~BeamformingFilter() { } -void mitk::BeamformingDMASFilter::GenerateInputRequestedRegion() +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::BeamformingDMASFilter::GenerateOutputInformation() +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()"); if (!m_Conf.partial) { m_Conf.bounds[0] = 0; m_Conf.bounds[1] = input->GetDimension(2) - 1; } - unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) - m_Conf.bounds[0] - ((input->GetDimension(2)-1) - m_Conf.bounds[1]) }; + unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) - m_Conf.bounds[0] - ((input->GetDimension(2) - 1) - m_Conf.bounds[1]) }; 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.RecordTime * 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::BeamformingDMASFilter::GenerateData() +void mitk::BeamformingFilter::GenerateData() { + GenerateOutputInformation(); mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; - double inputDim[2] = { input->GetDimension(0), input->GetDimension(1) / ((int)m_Conf.Photoacoustic + 1) }; + double inputDim[2] = { input->GetDimension(0), input->GetDimension(1) / ((int)m_Conf.Photoacoustic + 1) }; // if the photoacoustic option is used, we halve the image, as only the upper part of it contains any information double outputDim[2] = { output->GetDimension(0), output->GetDimension(1) }; const int apodArraySize = m_Conf.TransducerElements * 4; // set the resolution of the apodization array double* ApodWindow; // calculate the appropiate apodization window if (m_Conf.Apod == beamformingSettings::Apodization::Hann) { ApodWindow = VonHannFunction(apodArraySize); } else if (m_Conf.Apod == beamformingSettings::Apodization::Hamm) { ApodWindow = HammFunction(apodArraySize); } else { ApodWindow = BoxFunction(apodArraySize); } int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... for (int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i + m_Conf.bounds[0])); m_OutputData = new double[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; m_InputDataPuffer = new double[input->GetDimension(0)*input->GetDimension(1)]; // first, we convert any data to double, which we use by default if (input->GetPixelType().GetTypeAsString() == "scalar (double)" || input->GetPixelType().GetTypeAsString() == " (double)") { m_InputData = (double*)inputReadAccessor.GetData(); } else if (input->GetPixelType().GetTypeAsString() == "scalar (short)") { short* InputPuffer = (short*)inputReadAccessor.GetData(); for (int l = 0; l < inputDim[0]; ++l) { for (int s = 0; s < inputDim[1]; ++s) { m_InputDataPuffer[l*(short)inputDim[1] + s] = (double)InputPuffer[l*(short)inputDim[1] + s]; } } m_InputData = m_InputDataPuffer; } else if (input->GetPixelType().GetTypeAsString() == "scalar (float)") { float* InputPuffer = (float*)inputReadAccessor.GetData(); for (int l = 0; l < inputDim[0]; ++l) { for (int s = 0; s < inputDim[1]; ++s) { m_InputDataPuffer[l*(short)inputDim[1] + s] = (double)InputPuffer[l*(short)inputDim[1] + s]; } } m_InputData = m_InputDataPuffer; } else { MITK_INFO << "Could not determine pixel type"; return; } // 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; } } - const short max_threads = 12; - std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread - for (short line = 0; line < outputDim[0]; ++line) + if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { - threads[line] = std::thread(&BeamformingDMASFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); + for (short line = 0; line < outputDim[0]; ++line) + { + threads[line] = std::thread(&BeamformingFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); + } } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { - threads[line] = std::thread(&BeamformingDMASFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); + for (short line = 0; line < outputDim[0]; ++line) + { + threads[line] = std::thread(&BeamformingFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); + } } } - // wait for all lines to finish - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line].join(); - } - - /* - bool *threadfinished = new bool[max_threads]; - for (short r = 0; r < max_threads; ++r) - { - threadfinished[r] = false; - } - - short line = 0; - while(line < outputDim[0]) + else if(m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DMAS) { - //DMASSphericalLine(m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize) - //threads[line] = std::thread(&BeamformingDMASFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); - for (short n = 0; n < max_threads; ++n) + if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { - if (threadfinished[n]) + for (short line = 0; line < outputDim[0]; ++line) { - threads[n].join(); - threadfinished[n] = false; - //MITK_INFO << "thread " << n << " joined"; + threads[line] = std::thread(&BeamformingFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } - if (!threads[n].joinable()) + } + else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) + { + for (short line = 0; line < outputDim[0]; ++line) { - threads[n] = std::thread(&BeamformingDMASFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize, &threadfinished[n]); - ++line; - //MITK_INFO << "thread " << n << " created for line " << line - 1; - break; + threads[line] = std::thread(&BeamformingFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } } - for (short n = 0; n < max_threads; ++n) + // wait for all lines to finish + for (short line = 0; line < outputDim[0]; ++line) { - if (threads[n].joinable()) - { - threads[n].join(); - threadfinished[n] = false; - //MITK_INFO << "thread " << n << " joined"; - } + threads[line].join(); } - delete[] threadfinished;*/ - //threadpool... seems slower but I'll just keep it here for now output->SetSlice(m_OutputData, i); if (i % progInterval == 0) m_ProgressHandle((int)((i + 1) / (double)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; delete[] m_InputDataPuffer; - delete[] threads; m_OutputData = nullptr; m_InputData = nullptr; } // apply a bandpass filter, if requested if (m_Conf.UseBP) { m_ProgressHandle(100, "applying bandpass"); mitk::Image::Pointer BP = BandpassFilter(output); for (int i = 0; i < output->GetDimension(2); ++i) { mitk::ImageReadAccessor copy(BP, BP->GetSliceData(i)); output->SetSlice(copy.GetData(), i); } } m_TimeOfHeaderInitialization.Modified(); - delete[] ApodWindow; - auto end = std::chrono::high_resolution_clock::now(); - MITK_INFO << "DMAS Beamforming of " << output->GetDimension(2) << " Images completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; + MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } -void mitk::BeamformingDMASFilter::Configure(beamformingSettings settings) +void mitk::BeamformingFilter::Configure(beamformingSettings settings) { m_Conf = settings; } -mitk::Image::Pointer mitk::BeamformingDMASFilter::BandpassFilter(mitk::Image::Pointer data) +double* mitk::BeamformingFilter::VonHannFunction(int samples) +{ + double* ApodWindow = new double[samples]; + + for (int n = 0; n < samples; ++n) + { + ApodWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; + } + + return ApodWindow; +} + +double* mitk::BeamformingFilter::HammFunction(int samples) +{ + double* ApodWindow = new double[samples]; + + for (int n = 0; n < samples; ++n) + { + ApodWindow[n] = 0.54 - 0.46*cos(2 * M_PI*n / (samples - 1)); + } + + return ApodWindow; +} + +double* mitk::BeamformingFilter::BoxFunction(int samples) +{ + double* ApodWindow = new double[samples]; + + for (int n = 0; n < samples; ++n) + { + ApodWindow[n] = 1; + } + + return ApodWindow; +} + +void mitk::BeamformingFilter::DASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) +{ + double& inputS = inputDim[1]; + double& inputL = inputDim[0]; + + double& outputS = outputDim[1]; + double& outputL = outputDim[0]; + + short AddSample = 0; + short maxLine = 0; + short minLine = 0; + double delayMultiplicator = 0; + double l_i = 0; + double s_i = 0; + + double l = 0; + double x = 0; + double root = 0; + + double part = 0.07 * inputL; + double tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; + double apod_mult = 1; + + double 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; + + 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.0); + usedLines = (maxLine - minLine); + + apod_mult = apodArraySize / (maxLine - minLine); + + delayMultiplicator = pow((inputS / (m_Conf.RecordTime*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; + 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(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) +{ + double& inputS = inputDim[1]; + double& inputL = inputDim[0]; + + double& outputS = outputDim[1]; + double& outputL = outputDim[0]; + + short AddSample = 0; + short maxLine = 0; + short minLine = 0; + double delayMultiplicator = 0; + double l_i = 0; + double s_i = 0; + + double l = 0; + double x = 0; + double root = 0; + + double part = 0.07 * inputL; + double tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; + double apod_mult = 1; + + double mult = 0; + short usedLines = (maxLine - minLine); + + //exact delay + + l_i = (double)line / outputL * inputL; + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (double)sample / outputS * inputS; + + 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.0); + usedLines = (maxLine - minLine); + + apod_mult = apodArraySize / (maxLine - minLine); + + for (short l_s = minLine; l_s < maxLine; ++l_s) + { + AddSample = (int)sqrt( + pow(s_i, 2) + + + pow((inputS / (m_Conf.RecordTime*m_Conf.SpeedOfSound) * ((l_s - l_i)*m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) + ); + 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; + } +} + +mitk::Image::Pointer mitk::BeamformingFilter::BandpassFilter(mitk::Image::Pointer data) { // do a fourier transform, multiply with an appropriate window for the filter, and transform back typedef double 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 can not be applied after beamforming"; return data; } double singleVoxel = 1 / (m_Conf.RecordTime / data->GetDimension(1)) / 2 / 1000; double BoundHighPass = std::min(m_Conf.BPHighPass / singleVoxel, (double)data->GetDimension(1) / 2); double BoundLowPass = std::min(m_Conf.BPLowPass / singleVoxel, (double)data->GetDimension(1) / 2 - BoundHighPass); - int center1 = ((- BoundLowPass - BoundHighPass + data->GetDimension(1) / 2) / 2) + BoundLowPass; - int center2 = ((- BoundLowPass - BoundHighPass + data->GetDimension(1) / 2) / 2) + BoundHighPass + data->GetDimension(1) / 2; + int center1 = ((-BoundLowPass - BoundHighPass + data->GetDimension(1) / 2) / 2) + BoundLowPass; + int center2 = ((-BoundLowPass - BoundHighPass + data->GetDimension(1) / 2) / 2) + BoundHighPass + data->GetDimension(1) / 2; int width1 = -BoundLowPass - BoundHighPass + data->GetDimension(1) / 2; int width2 = -BoundLowPass - BoundHighPass + data->GetDimension(1) / 2; - /*MITK_INFO << "BHP " << BoundHighPass << " BLP " << BoundLowPass << "BPLP" << m_Conf.BPLowPass; MITK_INFO << "center1 " << center1 << " width1 " << width1; MITK_INFO << "center2 " << center2 << " width2 " << width2;*/ //DEBUG RealImageType::Pointer fftMultiplicator1 = BPFunction(data, width1, center1); RealImageType::Pointer fftMultiplicator2 = BPFunction(data, width2, center2); typedef itk::AddImageFilter AddImageFilterType; AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New(); addImageFilter->SetInput1(fftMultiplicator1); addImageFilter->SetInput2(fftMultiplicator2); typedef itk::FFTShiftImageFilter< RealImageType, RealImageType > FFTShiftFilterType; FFTShiftFilterType::Pointer fftShiftFilter = FFTShiftFilterType::New(); fftShiftFilter->SetInput(addImageFilter->GetOutput()); typedef itk::MultiplyImageFilter< ComplexImageType, RealImageType, ComplexImageType > MultiplyFilterType; MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); multiplyFilter->SetInput2(fftShiftFilter->GetOutput()); /*itk::ComplexToModulusImageFilter::Pointer toReal = itk::ComplexToModulusImageFilter::New(); toReal->SetInput(multiplyFilter->GetOutput()); return GrabItkImageMemory(addImageFilter->GetOutput());*/ //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::BeamformingDMASFilter::BPFunction(mitk::Image::Pointer reference, int width, int center) +itk::Image::Pointer mitk::BeamformingFilter::BPFunction(mitk::Image::Pointer reference, int width, int center) { - double* imageData = new double[reference->GetDimension(0)*reference->GetDimension(1)]; for (int sample = 0; sample < reference->GetDimension(1); ++sample) { imageData[reference->GetDimension(0)*sample] = 0; } /* // tukey window double alpha = m_Conf.BPFalloff; for (int n = 0; n < width; ++n) { - if (n <= (alpha*(width - 1)) / 2) - { - imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) - 1))) / 2; - } - else if (n >= (width - 1)*(1 - alpha / 2) && n <= (width - 1)) - { - imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; - } - else - { - imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = 1; - } + if (n <= (alpha*(width - 1)) / 2) + { + imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) - 1))) / 2; + } + else if (n >= (width - 1)*(1 - alpha / 2) && n <= (width - 1)) + { + imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; + } + else + { + imageData[reference->GetDimension(0)*(n + center - (int)(width / 2))] = 1; + } }*/ // Butterworth-Filter double d = center - width / 2; double l = center + width / 2; for (int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = 1 / (1 + pow(((double)center - (double)n) / ((double)width / 2), 2 * m_Conf.ButterworthOrder)); } // copy and paste to all lines for (int line = 1; line < reference->GetDimension(0); ++line) { for (int sample = 0; sample < reference->GetDimension(1); ++sample) { imageData[reference->GetDimension(0)*sample + line] = imageData[reference->GetDimension(0)*sample]; } } typedef itk::Image< double, 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; } -double* mitk::BeamformingDMASFilter::VonHannFunction(int samples) -{ - double* ApodWindow = new double[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; - } - - return ApodWindow; -} - -double* mitk::BeamformingDMASFilter::HammFunction(int samples) -{ - double* ApodWindow = new double[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = 0.54 - 0.46*cos(2 * M_PI*n / (samples - 1)); - } - - return ApodWindow; -} - -double* mitk::BeamformingDMASFilter::BoxFunction(int samples) -{ - double* ApodWindow = new double[samples]; - - for (int n = 0; n < samples; ++n) - { - ApodWindow[n] = 1; - } - - return ApodWindow; -} - -void mitk::BeamformingDMASFilter::DMASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) +void mitk::BeamformingFilter::DMASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; short AddSample1 = 0; short AddSample2 = 0; short maxLine = 0; short minLine = 0; double delayMultiplicator = 0; double l_i = 0; double s_i = 0; double l = 0; double x = 0; double root = 0; double part = 0.07 * inputL; double tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double apod_mult = 1; double 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; 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.0); usedLines = (maxLine - minLine); apod_mult = apodArraySize / (maxLine - minLine); delayMultiplicator = pow((inputS / (m_Conf.RecordTime*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); } for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < (short)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) { mult = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL] * apodisation[(short)((l_s2 - minLine)*apod_mult)] * input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL] * apodisation[(short)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = 10 * output[sample*(short)outputL + line] / (pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } -void mitk::BeamformingDMASFilter::DMASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) +void mitk::BeamformingFilter::DMASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; short AddSample1 = 0; short AddSample2 = 0; short maxLine = 0; short minLine = 0; double delayMultiplicator = 0; double l_i = 0; double s_i = 0; double l = 0; double x = 0; double root = 0; double part = 0.07 * inputL; double tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double apod_mult = 1; double mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; 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.0); usedLines = (maxLine - minLine); apod_mult = apodArraySize / (maxLine - minLine); //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((inputS / (m_Conf.RecordTime*m_Conf.SpeedOfSound) * ((minLine + l_s - l_i)*m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) ); } 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) { mult = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL] * apodisation[(int)((l_s2 - minLine)*apod_mult)] * input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL] * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = 10 * output[sample*(short)outputL + line] / (pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.h b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.h similarity index 78% rename from Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.h rename to Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.h index 02bdb9bd59..ddba028ece 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.h @@ -1,103 +1,109 @@ /*=================================================================== 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_DAS_FILTER -#define MITK_PHOTOACOUSTICS_BEAMFORMING_DAS_FILTER +#ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER +#define MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #include "mitkImageToImageFilter.h" #include namespace mitk { //##Documentation //## @brief //## @ingroup Process - class BeamformingDASFilter : public ImageToImageFilter + class BeamformingFilter : public ImageToImageFilter { public: - mitkClassMacro(BeamformingDASFilter, ImageToImageFilter); + mitkClassMacro(BeamformingFilter, ImageToImageFilter); itkFactorylessNewMacro(Self) itkCloneMacro(Self) struct beamformingSettings { double Pitch = 0.0003; // [m] float SpeedOfSound = 1540; // [m/s] unsigned int SamplesPerLine = 2048; unsigned int ReconstructionLines = 128; double RecordTime = 0.00006; // [s] unsigned int TransducerElements = 128; bool partial = false; unsigned int bounds[2] = { 0,0 }; enum DelayCalc {QuadApprox, Spherical}; DelayCalc DelayCalculationMethod = QuadApprox; enum Apodization {Hamm, Hann, Box}; Apodization Apod = Hann; + enum BeamformingAlgorithm {DMAS, DAS}; + BeamformingAlgorithm Algorithm = DAS; + double Angle = 10; bool Photoacoustic = true; double BPHighPass = 50; double BPLowPass = 50; bool UseBP = false; unsigned int ButterworthOrder = 8; }; void Configure(beamformingSettings settings); void SetProgressHandle(std::function progressHandle); protected: - BeamformingDASFilter(); + BeamformingFilter(); - ~BeamformingDASFilter(); + ~BeamformingFilter(); virtual void GenerateInputRequestedRegion() override; virtual void GenerateOutputInformation() override; virtual void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; std::function m_ProgressHandle; double* VonHannFunction(int samples); double* HammFunction(int samples); double* BoxFunction(int samples); void DASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize); void DASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize); + void DMASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize); + void DMASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const short& line, double* apodisation, const short& apodArraySize); + mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data); itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int width, int center); double* m_OutputData; double* m_InputData; double* m_InputDataPuffer; beamformingSettings m_Conf; }; } // namespace mitk -#endif MITK_PHOTOACOUSTICS_BEAMFORMING_DAS_FILTER +#endif MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt index 69b4fcb237..a1bd2a69cf 100644 --- a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt +++ b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt @@ -1,8 +1,8 @@ MITK_CREATE_MODULE( SUBPROJECTS - DEPENDS MitkCore MitkAlgorithmsExt + DEPENDS MitkCore MitkAlgorithmsExt MitkOpenCL #AUTOLOAD_WITH MitkCore INCLUDE_DIRS PUBLIC Algorithms/ITKUltrasound Algorithms Algorithms/OCL INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} PACKAGE_DEPENDS ITK|ITKFFT+ITKImageCompose+ITKImageIntensity ) diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index 69ecc38714..9096d9c260 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,6 +1,7 @@ set(CPP_FILES mitkPhotoacousticImage.cpp - Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp - Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp + Algorithms/mitkPhotoacousticBeamformingFilter.cpp + + Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp ) diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp index 993e2169eb..1ebc463364 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp @@ -1,287 +1,258 @@ /*=================================================================== 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 "itkBModeImageFilter.h" #include "itkPhotoacousticBModeImageFilter.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" -#include "mitkPhotoacousticBeamformingDASFilter.h" -#include "mitkPhotoacousticBeamformingDMASFilter.h" +#include "mitkPhotoacousticBeamformingFilter.h" #include #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" 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, bool UseLogFilter, float resampleSpacing) { // we use this seperate ApplyBmodeFilter Method for processing of two-dimensional images // 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::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; outputSize[0] = 256; outputSpacing[0] = itkImage->GetSpacing()[0] * (static_cast(inputSize[0]) / static_cast(outputSize[0])); outputSpacing[1] = resampleSpacing; outputSpacing[2] = 0.6; outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1]; typedef itk::IdentityTransform TransformType; resampleImageFilter->SetInput(bmode); resampleImageFilter->SetSize(outputSize); resampleImageFilter->SetOutputSpacing(outputSpacing); resampleImageFilter->SetTransform(TransformType::New()); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } /*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[3]) { typedef itk::Image< double, 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 inputSpacing = itkImage->GetSpacing(); itkFloatImageType::SpacingType outputSpacingItk; itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSizeItk = inputSizeItk; outputSizeItk[0] = outputSize[0]; outputSizeItk[1] = outputSize[1]; outputSizeItk[2] = outputSize[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] * (static_cast(inputSizeItk[2]) / static_cast(outputSizeItk[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) { double* inputData = new double[inputImage->GetDimension(0)*inputImage->GetDimension(1)*inputImage->GetDimension(2)]; double* outputData = new double[(inputImage->GetDimension(0) - left - right)*(inputImage->GetDimension(1) - above - below)*inputImage->GetDimension(2)]; 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) - above - below, inputImage->GetDimension(2) }; ImageReadAccessor acc(inputImage); // convert the data to double by default // as of now only those double, short, float are used at all... though it's easy to add other ones if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (double)") { inputData = (double*)acc.GetData(); } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (short)" || inputImage->GetPixelType().GetTypeAsString() == " (short)") { short* inputPuffer = (short*)acc.GetData(); for (int sl = 0; sl < inputDim[2]; ++sl) { for (int l = 0; l < inputDim[0]; ++l) { for (int s = 0; s < inputDim[1]; ++s) { inputData[l*inputDim[1] + s + sl*inputDim[0]*inputDim[1]] = (double)inputPuffer[l*inputDim[1] + s + sl*inputDim[0]*inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") { float* inputPuffer = (float*)acc.GetData(); for (int sl = 0; sl < inputDim[2]; ++sl) { for (int l = 0; l < inputDim[0]; ++l) { for (int s = 0; s < inputDim[1]; ++s) { inputData[l*inputDim[1] + s + sl*inputDim[0]*inputDim[1]] = (double)inputPuffer[l*inputDim[1] + s + sl*inputDim[0]*inputDim[1]]; } } } } else { MITK_INFO << "Could not determine pixel type"; } // copy the data into the cropped image for (int sl = 0; sl < outputDim[2]; ++sl) { for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] +sl*outputDim[0]*outputDim[1]] = inputData[(l + left) + (s + above)*(unsigned short)inputDim[0] + sl*inputDim[0]*inputDim[1]]; } } } mitk::Image::Pointer output = mitk::Image::New(); output->Initialize(mitk::MakeScalarPixelType(), 3, outputDim); output->SetSpacing(inputImage->GetGeometry()->GetSpacing()); for (int slice = 0; slice < outputDim[2]; ++slice) { output->SetSlice(&outputData[slice*outputDim[0] * outputDim[1]], slice); } return output; } -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamformingDAS(mitk::Image::Pointer inputImage, BeamformingDASFilter::beamformingSettings config, int cutoff, std::function progressHandle) +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingFilter::beamformingSettings config, int cutoff, std::function progressHandle) { // crop the image config.RecordTime = config.RecordTime - cutoff / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping progressHandle(0, "cropping image"); Image::Pointer croppedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0); Image::Pointer resizedImage = croppedImage; // resample the image in horizontal direction if (inputImage->GetDimension(0) != config.ReconstructionLines) { progressHandle(0, "resampling image"); auto begin = std::chrono::high_resolution_clock::now(); unsigned int dim[3] = { config.ReconstructionLines, croppedImage->GetDimension(1), croppedImage->GetDimension(2) }; resizedImage = ApplyResampling(croppedImage, dim); auto end = std::chrono::high_resolution_clock::now(); MITK_INFO << "Upsampling from " << inputImage->GetDimension(0) << " to " << config.ReconstructionLines << " lines completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } // perform the beamforming - BeamformingDASFilter::Pointer Beamformer = BeamformingDASFilter::New(); + BeamformingFilter::Pointer Beamformer = BeamformingFilter::New(); Beamformer->SetInput(resizedImage); Beamformer->Configure(config); Beamformer->SetProgressHandle(progressHandle); Beamformer->UpdateLargestPossibleRegion(); return Beamformer->GetOutput(); -} - -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamformingDMAS(mitk::Image::Pointer inputImage, BeamformingDMASFilter::beamformingSettings config, int cutoff, std::function progressHandle) -{ - // crop the image - config.RecordTime = config.RecordTime - cutoff / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping - progressHandle(0, "cropping image"); - Image::Pointer croppedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0); - Image::Pointer resizedImage = croppedImage; - - // resample the image in horizontal direction - if(inputImage->GetDimension(0) != config.ReconstructionLines) - { - progressHandle(0, "resampling image"); - auto begin = std::chrono::high_resolution_clock::now(); - unsigned int dim[3] = { config.ReconstructionLines, croppedImage->GetDimension(1), croppedImage->GetDimension(2) }; - resizedImage = ApplyResampling(croppedImage, dim); - auto end = std::chrono::high_resolution_clock::now(); - MITK_INFO << "Upsampling from " << inputImage->GetDimension(0) << " to " << config.ReconstructionLines << " lines completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; - } - - // perform the beamforming - BeamformingDMASFilter::Pointer Beamformer = BeamformingDMASFilter::New(); - Beamformer->SetInput(resizedImage); - Beamformer->Configure(config); - Beamformer->SetProgressHandle(progressHandle); - Beamformer->UpdateLargestPossibleRegion(); - return Beamformer->GetOutput(); -} +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h index 1f931c35ee..eea9ba1ca3 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h @@ -1,51 +1,49 @@ /*=================================================================== 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 "Algorithms/mitkPhotoacousticBeamformingDASFilter.h" -#include "Algorithms/mitkPhotoacousticBeamformingDMASFilter.h" +#include "mitkPhotoacousticBeamformingFilter.h" #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticImage : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticImage, itk::Object); itkFactorylessNewMacro(Self); mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, 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[3]); - mitk::Image::Pointer ApplyBeamformingDAS(mitk::Image::Pointer inputImage, BeamformingDASFilter::beamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); - mitk::Image::Pointer ApplyBeamformingDMAS(mitk::Image::Pointer inputImage, BeamformingDMASFilter::beamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); + mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingFilter::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); protected: PhotoacousticImage(); virtual ~PhotoacousticImage(); }; } // namespace mitk #endif /* mitkPhotoacousticImage_H_HEADER_INCLUDED */ 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 3e1eb9f930..0b2e61b22d 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,683 +1,639 @@ /*=================================================================== 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 //mitk image #include #include "mitkPhotoacousticImage.h" -#include "mitkPhotoacousticBeamformingDASFilter.h" -#include "mitkPhotoacousticBeamformingDMASFilter.h" +#include "mitkPhotoacousticBeamformingFilter.h" //other #include #include const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false) { 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.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UseImageSpacing())); connect(m_Controls.ScanDepth, SIGNAL(valueChanged(double)), this, SLOT(UpdateFrequency())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(UpdateFrequency())); connect(m_Controls.Samples, SIGNAL(valueChanged(int)), this, SLOT(UpdateFrequency())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UpdateFrequency())); connect(m_Controls.UseBP, SIGNAL(clicked()), this, SLOT(UseBandpass())); connect(m_Controls.boundLow, SIGNAL(valueChanged(int)), this, SLOT(UpdateBounds())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpdateBounds())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(UpdateBounds())); 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); UseImageSpacing(); UseBandpass(); } 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::finished, thread, &QObject::deleteLater); - thread->setBeamformingAlgorithm(m_CurrentBeamformingAlgorithm); - thread->setConfigs(DMASconfig, DASconfig); + thread->setConfig(BFconfig); thread->setInputImage(image); 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 (m_CurrentBeamformingAlgorithm == BeamformingAlgorithms::DAS) + if (BFconfig.Algorithm == mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS) newNodeName << "DAS bf, "; - else if (m_CurrentBeamformingAlgorithm == BeamformingAlgorithms::DMAS) + else if (BFconfig.Algorithm == mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DMAS) newNodeName << "DMAS bf, "; - if (DASconfig.DelayCalculationMethod == mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::QuadApprox) + if (BFconfig.DelayCalculationMethod == mitk::BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox) newNodeName << "q. delay"; - if (DASconfig.DelayCalculationMethod == mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::Spherical) + if (BFconfig.DelayCalculationMethod == mitk::BeamformingFilter::beamformingSettings::DelayCalc::Spherical) newNodeName << "s. delay"; 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.ProgressInfo->setVisible(false); m_Controls.buttonApplyBeamforming->setText("Apply Beamforming"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->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); thread->setConfig(m_UseLogfilter, m_ResampleSpacing); thread->setInputImage(image); 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::UpdateBounds() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); - DASconfig.partial = false; - DMASconfig.partial = false; - MITK_INFO << "ASFASF"; + BFconfig.partial = false; return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); - DASconfig.partial = true; - DMASconfig.partial = true; + BFconfig.partial = true; } if(m_Controls.boundLow->value()>m_Controls.boundHigh->value()) { MITK_INFO << "high bound < low bound -> setting both to beamform only first slice"; m_Controls.boundLow->setValue(0); m_Controls.boundHigh->setValue(0); - DASconfig.bounds[0] = 0; - DMASconfig.bounds[0] = 0; - DASconfig.bounds[1] = 0; - DMASconfig.bounds[1] = 0; + BFconfig.bounds[0] = 0; + BFconfig.bounds[1] = 0; } else { - DASconfig.bounds[0] = m_Controls.boundLow->value(); - DMASconfig.bounds[0] = m_Controls.boundLow->value(); - DASconfig.bounds[1] = m_Controls.boundHigh->value(); - DMASconfig.bounds[1] = m_Controls.boundHigh->value(); + BFconfig.bounds[0] = m_Controls.boundLow->value(); + BFconfig.bounds[1] = m_Controls.boundHigh->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::UpdateFrequency() { - DMASconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] - DMASconfig.SamplesPerLine = m_Controls.Samples->value(); + BFconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] + BFconfig.SamplesPerLine = m_Controls.Samples->value(); QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected if (!m_Controls.UseImageSpacing->isChecked()) UpdateRecordTime(mitk::Image::New()); 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(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); } if (image) UpdateRecordTime(image); } std::stringstream frequency; - frequency << 1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000 / 1000000; //[MHz] + frequency << 1 / (BFconfig.RecordTime / BFconfig.SamplesPerLine) * BFconfig.SamplesPerLine / 2 / 2 / 1000 / 1000000; //[MHz] frequency << "MHz"; 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 ); UpdateFrequency(); return; } } m_Controls.labelWarning->setVisible( true ); m_Controls.buttonApplyBModeFilter->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()) - m_CurrentBeamformingAlgorithm = BeamformingAlgorithms::DAS; + BFconfig.Algorithm = mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) - m_CurrentBeamformingAlgorithm = BeamformingAlgorithms::DMAS; + BFconfig.Algorithm = mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DMAS; if ("Quad. Approx." == m_Controls.DelayCalculation->currentText()) { - DASconfig.DelayCalculationMethod = mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::QuadApprox; - DMASconfig.DelayCalculationMethod = mitk::BeamformingDMASFilter::beamformingSettings::DelayCalc::QuadApprox; + BFconfig.DelayCalculationMethod = mitk::BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Wave" == m_Controls.DelayCalculation->currentText()) { - DASconfig.DelayCalculationMethod = mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::Spherical; - DMASconfig.DelayCalculationMethod = mitk::BeamformingDMASFilter::beamformingSettings::DelayCalc::Spherical; + BFconfig.DelayCalculationMethod = mitk::BeamformingFilter::beamformingSettings::DelayCalc::Spherical; } if ("Von Hann" == m_Controls.Apodization->currentText()) { - DASconfig.Apod = mitk::BeamformingDASFilter::beamformingSettings::Apodization::Hann; - DMASconfig.Apod = mitk::BeamformingDMASFilter::beamformingSettings::Apodization::Hann; + BFconfig.Apod = mitk::BeamformingFilter::beamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { - DASconfig.Apod = mitk::BeamformingDASFilter::beamformingSettings::Apodization::Hamm; - DMASconfig.Apod = mitk::BeamformingDMASFilter::beamformingSettings::Apodization::Hamm; + BFconfig.Apod = mitk::BeamformingFilter::beamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { - DASconfig.Apod = mitk::BeamformingDASFilter::beamformingSettings::Apodization::Box; - DMASconfig.Apod = mitk::BeamformingDMASFilter::beamformingSettings::Apodization::Box; + BFconfig.Apod = mitk::BeamformingFilter::beamformingSettings::Apodization::Box; } - DASconfig.Pitch = m_Controls.Pitch->value() / 1000; // [m] - DASconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] - DASconfig.SamplesPerLine = m_Controls.Samples->value(); - DASconfig.ReconstructionLines = m_Controls.Lines->value(); - DASconfig.TransducerElements = m_Controls.ElementCount->value(); - DASconfig.Angle = m_Controls.Angle->value(); // [deg] - DASconfig.BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] - DASconfig.BPLowPass = 1000000 * (1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000 / 1000000 - m_Controls.BPlow->value()); // [Hz] - DASconfig.ButterworthOrder = m_Controls.BPFalloff->value(); - DASconfig.UseBP = m_Controls.UseBP->isChecked(); - - DMASconfig.Pitch = m_Controls.Pitch->value() / 1000; // [m] - DMASconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] - DMASconfig.SamplesPerLine = m_Controls.Samples->value(); - DMASconfig.ReconstructionLines = m_Controls.Lines->value(); - DMASconfig.TransducerElements = m_Controls.ElementCount->value(); - DMASconfig.Angle = m_Controls.Angle->value(); //[deg] - DMASconfig.BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] - DMASconfig.BPLowPass = 1000000 * (1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000 /1000000 - m_Controls.BPlow->value()); // [Hz] - DMASconfig.ButterworthOrder = m_Controls.BPFalloff->value(); - DMASconfig.UseBP = m_Controls.UseBP->isChecked(); + 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.Angle = m_Controls.Angle->value(); // [deg] + BFconfig.BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] + BFconfig.BPLowPass = 1000000 * (1 / (BFconfig.RecordTime / BFconfig.SamplesPerLine) * BFconfig.SamplesPerLine / 2 / 2 / 1000 / 1000000 - m_Controls.BPlow->value()); // [Hz] + BFconfig.ButterworthOrder = m_Controls.BPFalloff->value(); + BFconfig.UseBP = m_Controls.UseBP->isChecked(); UpdateRecordTime(image); UpdateBounds(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image - double maxFrequency = 1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000; // [Hz] + double maxFrequency = 1 / (BFconfig.RecordTime / BFconfig.SamplesPerLine) * BFconfig.SamplesPerLine / 2 / 2 / 1000; // [Hz] - if (DMASconfig.BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) + if (BFconfig.BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); - DMASconfig.BPLowPass = 0; - DASconfig.BPLowPass = 0; + BFconfig.BPLowPass = 0; } - if (DMASconfig.BPLowPass < 0 && m_Controls.UseBP->isChecked()) + if (BFconfig.BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); - DMASconfig.BPLowPass = 0; - DASconfig.BPLowPass = 0; + BFconfig.BPLowPass = 0; } - if (DMASconfig.BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) + if (BFconfig.BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); - DMASconfig.BPHighPass = 0; - DASconfig.BPHighPass = 0; + BFconfig.BPHighPass = 0; } - if (DMASconfig.BPHighPass > maxFrequency - DMASconfig.BPLowPass) + if (BFconfig.BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); - DMASconfig.BPHighPass = 0; - DASconfig.BPHighPass = 0; - DMASconfig.BPLowPass = 0; - DASconfig.BPLowPass = 0; + BFconfig.BPHighPass = 0; + BFconfig.BPLowPass = 0; } } void PAImageProcessing::UpdateRecordTime(mitk::Image::Pointer image) { if (m_Controls.UseImageSpacing->isChecked()) { - DASconfig.RecordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000; // [s] - DMASconfig.RecordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000; // [s] - MITK_INFO << "Calculated Scan Depth of " << DASconfig.RecordTime * DASconfig.SpeedOfSound * 100 << "cm"; + BFconfig.RecordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000; // [s] + MITK_INFO << "Calculated Scan Depth of " << BFconfig.RecordTime * BFconfig.SpeedOfSound * 100 << "cm"; } else { - DASconfig.RecordTime = m_Controls.ScanDepth->value() / 1000 / DASconfig.SpeedOfSound; // [s] - DMASconfig.RecordTime = m_Controls.ScanDepth->value() / 1000 / DMASconfig.SpeedOfSound; // [s] + BFconfig.RecordTime = m_Controls.ScanDepth->value() / 1000 / BFconfig.SpeedOfSound; // [s] } if ("US Image" == m_Controls.ImageType->currentText()) { if (m_Controls.UseImageSpacing->isChecked()) { - DASconfig.RecordTime = DASconfig.RecordTime / 2; // [s] - DMASconfig.RecordTime = DMASconfig.RecordTime / 2; // [s] + BFconfig.RecordTime = BFconfig.RecordTime / 2; // [s] } - DASconfig.Photoacoustic = false; - DMASconfig.Photoacoustic = false; + BFconfig.Photoacoustic = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { - DASconfig.Photoacoustic = true; - DMASconfig.Photoacoustic = true; + BFconfig.Photoacoustic = true; } } void PAImageProcessing::EnableControls() { m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.Cutoff->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); UseBandpass(); 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.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.Cutoff->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); 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); } UpdateFrequency(); } void PAImageProcessing::UseBandpass() { if (m_Controls.UseBP->isChecked()) { m_Controls.BPhigh->setEnabled(true); m_Controls.BPlow->setEnabled(true); m_Controls.BPFalloff->setEnabled(true); } else { m_Controls.BPhigh->setDisabled(true); m_Controls.BPlow->setDisabled(true); m_Controls.BPFalloff->setDisabled(true); } UpdateFrequency(); } void BeamformingThread::run() { mitk::Image::Pointer resultImage; mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; - if (m_CurrentBeamformingAlgorithm == PAImageProcessing::BeamformingAlgorithms::DAS) - resultImage = filterbank->ApplyBeamformingDAS(m_InputImage, m_DASconfig, m_Cutoff, progressHandle); - else if (m_CurrentBeamformingAlgorithm == PAImageProcessing::BeamformingAlgorithms::DMAS) - resultImage = filterbank->ApplyBeamformingDMAS(m_InputImage, m_DMASconfig, m_Cutoff, progressHandle); + resultImage = filterbank->ApplyBeamforming(m_InputImage, m_BFconfig, m_Cutoff, progressHandle); emit result(resultImage); } -void BeamformingThread::setConfigs(mitk::BeamformingDMASFilter::beamformingSettings DMASconfig, mitk::BeamformingDASFilter::beamformingSettings DASconfig) +void BeamformingThread::setConfig(mitk::BeamformingFilter::beamformingSettings BFconfig) { - m_DMASconfig = DMASconfig; - m_DASconfig = DASconfig; -} - -void BeamformingThread::setBeamformingAlgorithm(PAImageProcessing::BeamformingAlgorithms beamformingAlgorithm) -{ - m_CurrentBeamformingAlgorithm = beamformingAlgorithm; + m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BeamformingThread::setCutoff(int cutoff) { m_Cutoff = cutoff; } void BmodeThread::run() { mitk::Image::Pointer resultImage; mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); resultImage = filterbank->ApplyBmodeFilter(m_InputImage, m_UseLogfilter, m_ResampleSpacing); emit result(resultImage); } void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing) { m_UseLogfilter = useLogfilter; m_ResampleSpacing = resampleSpacing; } void BmodeThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } \ No newline at end of file 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 b79e3f320c..d09a380062 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,139 +1,129 @@ /*=================================================================== 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 "ui_PAImageProcessingControls.h" -#include "mitkPhotoacousticBeamformingDASFilter.h" -#include "mitkPhotoacousticBeamformingDMASFilter.h" +#include "mitkPhotoacousticBeamformingFilter.h" Q_DECLARE_METATYPE(mitk::Image::Pointer) Q_DECLARE_METATYPE(std::string) 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(); - enum BeamformingAlgorithms { DAS, DMAS }; - protected slots: /// \brief Called when the user clicks the GUI button void UpdateBounds(); void UseResampling(); void UseLogfilter(); void SetResampling(); void UseImageSpacing(); void UpdateFrequency(); void UpdateRecordTime(mitk::Image::Pointer image); void UseBandpass(); void HandleBeamformingResults(mitk::Image::Pointer image); void StartBeamformingThread(); void HandleBmodeResults(mitk::Image::Pointer image); void StartBmodeThread(); void UpdateProgress(int progress, std::string progressInfo); protected: virtual void CreateQtPartControl(QWidget *parent) override; virtual void SetFocus() override; /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( berry::IWorkbenchPart::Pointer source, const QList& nodes ) override; Ui::PAImageProcessingControls m_Controls; float m_ResampleSpacing; bool m_UseLogfilter; std::string m_OldNodeName; - mitk::BeamformingDMASFilter::beamformingSettings DMASconfig; - mitk::BeamformingDASFilter::beamformingSettings DASconfig; - - BeamformingAlgorithms m_CurrentBeamformingAlgorithm; + mitk::BeamformingFilter::beamformingSettings BFconfig; void UpdateBFSettings(mitk::Image::Pointer image); void EnableControls(); void DisableControls(); }; class BeamformingThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); void updateProgress(int, std::string); public: - void setConfigs(mitk::BeamformingDMASFilter::beamformingSettings DMASconfig, mitk::BeamformingDASFilter::beamformingSettings DASconfig); - void setBeamformingAlgorithm(PAImageProcessing::BeamformingAlgorithms beamformingAlgorithm); + void setConfig(mitk::BeamformingFilter::beamformingSettings BFconfig); void setInputImage(mitk::Image::Pointer image); void setCutoff(int cutoff); protected: - mitk::BeamformingDMASFilter::beamformingSettings m_DMASconfig; - mitk::BeamformingDASFilter::beamformingSettings m_DASconfig; - - PAImageProcessing::BeamformingAlgorithms m_CurrentBeamformingAlgorithm; + mitk::BeamformingFilter::beamformingSettings m_BFconfig; mitk::Image::Pointer m_InputImage; int m_Cutoff; }; class BmodeThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(bool useLogfilter, double resampleSpacing); void setInputImage(mitk::Image::Pointer image); protected: mitk::Image::Pointer m_InputImage; bool m_UseLogfilter; double m_ResampleSpacing; }; #endif // PAImageProcessing_h