diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp index 0044adbcf5..0d6b055f0d 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp @@ -1,180 +1,323 @@ /*=================================================================== 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("DASQuadratic.cl"); this->AddSourceFile("DMASQuadratic.cl"); this->AddSourceFile("DASspherical.cl"); this->AddSourceFile("DMASspherical.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, m_OutputDim); } catch (const mitk::Exception& e) { MITK_ERROR << "Catched exception while initializing filter: " << e.what(); return; } if (m_Apodisation == nullptr) { MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; m_Apodisation = new float[1]; m_Apodisation[0] = 1; m_ApodArraySize = 1; } us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); cl_mem cl_input = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, m_Apodisation, &clErr); CHECK_OCL_ERR(clErr); // set kernel arguments clErr = clSetKernelArg( this->m_PixelCalculation, 2, sizeof(cl_mem), &cl_input); clErr |= clSetKernelArg( this->m_PixelCalculation, 3, sizeof(cl_ushort), &(this->m_ApodArraySize) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 4, sizeof(cl_ushort), &(this->m_OutputDim[1]) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_OutputDim[0]) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 6, sizeof(cl_float), &(this->m_SpeedOfSound) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_RecordTime) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 8, sizeof(cl_float), &(this->m_Pitch) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 9, sizeof(cl_float), &(this->m_Angle) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 10, sizeof(cl_ushort), &(this->m_PAImage)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_ushort), &(this->m_TransducerElements)); CHECK_OCL_ERR( clErr ); // execute the filter on a 3D NDRange this->ExecuteKernel( m_PixelCalculation, 3); // signalize the GPU-side data changed m_Output->Modified( GPU_DATA ); } us::Module *mitk::PhotoacousticOCLBeamformer::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBeamformer::Initialize() { bool buildErr = true; cl_int clErr = 0; if ( OclFilter::Initialize() ) { switch (m_Algorithm) { case BeamformingAlgorithm::DASQuad: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASQuad", &clErr); break; } case BeamformingAlgorithm::DMASQuad: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMASQuad", &clErr); break; } case BeamformingAlgorithm::DASSphe: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASSphe", &clErr); break; } case BeamformingAlgorithm::DMASSphe: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMASSphe", &clErr); break; } } 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); - m_InputS = image->GetDimension(1); - m_InputL = image->GetDimension(0); } mitk::Image::Pointer mitk::PhotoacousticOCLBeamformer::GetOutput() { if (m_Output->IsModified(GPU_DATA)) { - void* pData = m_Output->TransferDataToCPU(m_CommandQue); + if (false)//(m_Algorithm == BeamformingAlgorithm::DMASQuad || m_Algorithm == BeamformingAlgorithm::DMASSphe) + { + mitk::PhotoacousticPixelSum::Pointer m_oclFilter = mitk::PhotoacousticPixelSum::New(); + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + resources->GetContext(); + + const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->GetSlicedGeometry(); + const unsigned int dimension = 3; + unsigned int* dimensions = m_OutputDim; + void* pData = m_Output->TransferDataToCPU(m_CommandQue); + + mitk::Image::Pointer temp = mitk::Image::New(); + temp->Initialize(this->GetOutputType(), dimension, dimensions); + temp->SetSpacing(p_slg->GetSpacing()); + temp->SetGeometry(m_Input->GetMITKImage()->GetGeometry()); + + m_oclFilter->SetInput(temp); + m_oclFilter->SetOutputDim(m_OutputDim); + m_oclFilter->SetData(pData); + m_oclFilter->SetImagesToAdd(m_OutputDim[0]); + m_oclFilter->Update(); + + void* out = m_oclFilter->GetRawOutput(); - const unsigned int dimension = 3; - unsigned int* dimensions = m_OutputDim; + MITK_DEBUG << "Creating new MITK Image."; - const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->GetSlicedGeometry(); + m_Output->GetMITKImage()->Initialize(this->GetOutputType(), dimension, dimensions); + m_Output->GetMITKImage()->SetSpacing(p_slg->GetSpacing()); + m_Output->GetMITKImage()->SetGeometry(m_Input->GetMITKImage()->GetGeometry()); + m_Output->GetMITKImage()->SetImportVolume(out, 0, 0, mitk::Image::ReferenceMemory); + } + else + { + void* pData = m_Output->TransferDataToCPU(m_CommandQue); + + const unsigned int dimension = 3; + unsigned int* dimensions = m_OutputDim; - MITK_DEBUG << "Creating new MITK Image."; + const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->GetSlicedGeometry(); - m_Output->GetMITKImage()->Initialize(this->GetOutputType(), dimension, dimensions); - m_Output->GetMITKImage()->SetSpacing(p_slg->GetSpacing()); - m_Output->GetMITKImage()->SetGeometry(m_Input->GetMITKImage()->GetGeometry()); - m_Output->GetMITKImage()->SetImportVolume(pData, 0, 0, mitk::Image::ReferenceMemory); + MITK_DEBUG << "Creating new MITK Image."; + + m_Output->GetMITKImage()->Initialize(this->GetOutputType(), dimension, dimensions); + m_Output->GetMITKImage()->SetSpacing(p_slg->GetSpacing()); + m_Output->GetMITKImage()->SetGeometry(m_Input->GetMITKImage()->GetGeometry()); + m_Output->GetMITKImage()->SetImportVolume(pData, 0, 0, mitk::Image::ReferenceMemory); + } } MITK_DEBUG << "Image Initialized."; return m_Output->GetMITKImage(); +} + +mitk::PhotoacousticPixelSum::PhotoacousticPixelSum() + : m_PixelCalculation(NULL) +{ + this->AddSourceFile("PixelSum.cl"); + + this->m_FilterID = "PixelSum"; +} + +mitk::PhotoacousticPixelSum::~PhotoacousticPixelSum() +{ + if (this->m_PixelCalculation) + { + clReleaseKernel(m_PixelCalculation); + } +} + +void mitk::PhotoacousticPixelSum::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::PhotoacousticPixelSum::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; + } + + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + cl_context gpuContext = resources->GetContext(); + + unsigned int size = m_OutputDim[0] * m_OutputDim[1] * m_OutputDim[2] * m_Images; + + cl_mem cl_input = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * size, m_Data, &clErr); + CHECK_OCL_ERR(clErr); + + // set kernel arguments + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &cl_input); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_ushort), &(this->m_Images)); + + CHECK_OCL_ERR(clErr); + + // execute the filter on a 3D NDRange + this->ExecuteKernel(m_PixelCalculation, 3); + + // signalize the GPU-side data changed + m_Output->Modified(GPU_DATA); +} + +us::Module *mitk::PhotoacousticPixelSum::GetModule() +{ + return us::GetModuleContext()->GetModule(); +} + +bool mitk::PhotoacousticPixelSum::Initialize() +{ + bool buildErr = true; + cl_int clErr = 0; + + if (OclFilter::Initialize()) + { + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckSum", &clErr); + buildErr |= CHECK_OCL_ERR(clErr); + } + return (OclFilter::IsInitialized() && buildErr); +} + +void mitk::PhotoacousticPixelSum::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); +} + +void* mitk::PhotoacousticPixelSum::GetRawOutput() +{ + if (m_Output->IsModified(GPU_DATA)) + { + void* pData = m_Output->TransferDataToCPU(m_CommandQue); + return pData; + } + else + MITK_ERROR << "No Images Summed"; } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h index ab9d44c805..55f31ab4fa 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h @@ -1,128 +1,198 @@ /*=================================================================== 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 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); mitk::Image::Pointer GetOutput(); /** 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]; } void SetApodisation(float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } enum BeamformingAlgorithm { DASQuad, DMASQuad, DASSphe, DMASSphe }; void SetAlgorithm(BeamformingAlgorithm algorithm, bool PA) { m_Algorithm = algorithm; m_PAImage = PA; } - void SetBeamformingParameters(float SpeedOfSound, float RecordTime, float Pitch, float Angle, bool PAImage) + void SetBeamformingParameters(float SpeedOfSound, float RecordTime, float Pitch, float Angle, bool PAImage, unsigned short transducerElements) { m_SpeedOfSound = SpeedOfSound; m_RecordTime = RecordTime; m_Pitch = Pitch; m_Angle = Angle; m_PAImage = PAImage; + m_TransducerElements = transducerElements; } 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]; - unsigned short m_InputS; - unsigned short m_InputL; float* m_Apodisation; unsigned short m_ApodArraySize; BeamformingAlgorithm m_Algorithm; unsigned short m_PAImage; float m_SpeedOfSound; float m_RecordTime; float m_Pitch; float m_Angle; + unsigned short m_TransducerElements; +}; + +class PhotoacousticPixelSum : public OclImageToImageFilter, public itk::Object +{ + +public: + mitkClassMacroItkParent(PhotoacousticPixelSum, 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); + + void* GetRawOutput(); + + /** 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]; + } + + void SetImagesToAdd(unsigned short number) + { + m_Images = number; + } + + void SetData(void* data) + { + m_Data = (float*)data; + } + +protected: + + /** Constructor */ + PhotoacousticPixelSum(); + + /** Destructor */ + virtual ~PhotoacousticPixelSum(); + + /** 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]; + float* m_Data; + unsigned short m_Images; }; } #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp index d9eb0ccaef..ceca2868b7 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp @@ -1,755 +1,756 @@ /*=================================================================== 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 "mitkPhotoacousticBeamformingFilter.h" #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include #include "itkFFT1DComplexConjugateToRealImageFilter.h" #include "itkFFT1DRealToComplexConjugateImageFilter.h" #include "mitkImageCast.h" #include // needed itk image filters #include "mitkITKImageImport.h" #include "itkFFTShiftImageFilter.h" #include "itkMultiplyImageFilter.h" #include "itkComplexToModulusImageFilter.h" #include mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); //GenerateTimeInInputRegion(output, input); } void mitk::BeamformingFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; itkDebugMacro(<< "GenerateOutputInformation()"); unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf.Pitch * m_Conf.TransducerElements * 1000 / m_Conf.ReconstructionLines; spacing[1] = m_Conf.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::BeamformingFilter::GenerateData() { GenerateOutputInformation(); mitk::Image::Pointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; float 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 float outputDim[2] = { output->GetDimension(0), output->GetDimension(1) }; unsigned int oclOutputDim[3] = { output->GetDimension(0), output->GetDimension(1), output->GetDimension(2) }; const int apodArraySize = m_Conf.TransducerElements * 4; // set the resolution of the apodization array float* 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_OutputData = new float[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; m_InputDataPuffer = new float[input->GetDimension(0)*input->GetDimension(1)]; // first, we convert any data to float, which we use by default if (input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)") { m_InputData = (float*)inputReadAccessor.GetData(); } else { MITK_INFO << "Pixel type is not float, abort"; return; } // fill the image with zeros for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } } else if(m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, 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) / (float)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; delete[] m_InputDataPuffer; m_OutputData = nullptr; m_InputData = nullptr; } */ mitk::PhotoacousticOCLBeamformer::Pointer m_oclFilter = mitk::PhotoacousticOCLBeamformer::New(); us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); resources->GetContext(); try { + MITK_INFO << "afas" << input->GetPixelType().GetTypeAsString(); m_oclFilter->SetInput(input); if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DAS) { if(m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DASQuad, true); else if(m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DASSphe, true); } else if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DMASQuad, true); else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DMASSphe, true); } m_oclFilter->SetApodisation(ApodWindow, apodArraySize); m_oclFilter->SetOutputDim(oclOutputDim); - m_oclFilter->SetBeamformingParameters(m_Conf.SpeedOfSound, m_Conf.RecordTime, m_Conf.Pitch, m_Conf.Angle, m_Conf.Photoacoustic); + m_oclFilter->SetBeamformingParameters(m_Conf.SpeedOfSound, m_Conf.RecordTime, m_Conf.Pitch, m_Conf.Angle, m_Conf.Photoacoustic, m_Conf.TransducerElements); m_oclFilter->Update(); auto out = m_oclFilter->GetOutput(); for (int i = 0; i < output->GetDimension(2); ++i) { mitk::ImageReadAccessor copy(out, out->GetSliceData(i)); output->SetSlice(copy.GetData(), i); } } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; } // 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 << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } void mitk::BeamformingFilter::Configure(beamformingSettings settings) { m_Conf = settings; } float* mitk::BeamformingFilter::VonHannFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; } return ApodWindow; } float* mitk::BeamformingFilter::HammFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 0.54 - 0.46*cos(2 * M_PI*n / (samples - 1)); } return ApodWindow; } float* mitk::BeamformingFilter::BoxFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 1; } return ApodWindow; } void mitk::BeamformingFilter::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float l = 0; float x = 0; float root = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = 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(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = 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 float PixelType; typedef itk::Image< PixelType, 3 > RealImageType; RealImageType::Pointer image; mitk::CastToItkImage(data, image); typedef itk::FFT1DRealToComplexConjugateImageFilter ForwardFFTFilterType; typedef ForwardFFTFilterType::OutputImageType ComplexImageType; ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New(); forwardFFTFilter->SetInput(image); forwardFFTFilter->SetDirection(1); try { forwardFFTFilter->UpdateOutputInformation(); } catch (itk::ExceptionObject & error) { std::cerr << "Error: " << error << std::endl; MITK_WARN << "Bandpass can not be applied after beamforming"; return data; } float singleVoxel = 1 / (m_Conf.RecordTime / data->GetDimension(1)) / 2 / 1000; float BoundHighPass = std::min(m_Conf.BPHighPass / singleVoxel, (float)data->GetDimension(1) / 2); float BoundLowPass = std::min(m_Conf.BPLowPass / singleVoxel, (float)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::BeamformingFilter::BPFunction(mitk::Image::Pointer reference, int width, int center) { float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; for (int sample = 0; sample < reference->GetDimension(1); ++sample) { imageData[reference->GetDimension(0)*sample] = 0; } /* // tukey window float 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 float d = center - width / 2; float l = center + width / 2; for (int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = 1 / (1 + pow(((float)center - (float)n) / ((float)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< float, 3U > ImageType; ImageType::RegionType region; ImageType::IndexType start; start.Fill(0); region.SetIndex(start); ImageType::SizeType size; size[0] = reference->GetDimension(0); size[1] = reference->GetDimension(1); size[2] = reference->GetDimension(2); region.SetSize(size); ImageType::SpacingType SpacingItk; SpacingItk[0] = reference->GetGeometry()->GetSpacing()[0]; SpacingItk[1] = reference->GetGeometry()->GetSpacing()[1]; SpacingItk[2] = reference->GetGeometry()->GetSpacing()[2]; ImageType::Pointer image = ImageType::New(); image->SetRegions(region); image->Allocate(); image->FillBuffer(itk::NumericTraits::Zero); image->SetSpacing(SpacingItk); ImageType::IndexType pixelIndex; for (ImageType::IndexValueType slice = 0; slice < reference->GetDimension(2); ++slice) { for (ImageType::IndexValueType line = 0; line < reference->GetDimension(0); ++line) { for (ImageType::IndexValueType sample = 0; sample < reference->GetDimension(1); ++sample) { pixelIndex[0] = line; pixelIndex[1] = sample; pixelIndex[2] = slice; image->SetPixel(pixelIndex, imageData[line + sample*reference->GetDimension(0)]); } } } delete[] imageData; return image; } void mitk::BeamformingFilter::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample1 = 0; short AddSample2 = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float l = 0; float x = 0; float root = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = 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::BeamformingFilter::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample1 = 0; short AddSample2 = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float l = 0; float x = 0; float root = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; float apod_mult = 1; float 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.0f); 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/Resources/DASQuadratic.cl b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl index 2918efcf8a..a7f6226a3a 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl @@ -1,77 +1,78 @@ /*=================================================================== 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 ckDASQuad( __read_only image3d_t dSource, // input image __global float* dDest, // output buffer __global float* apodArray, unsigned short apodArraySize, unsigned short outputS, unsigned short outputL, float SpeedOfSound, float RecordTime, float Pitch, float Angle, - unsigned short PAImage // parameters + unsigned short PAImage, + unsigned short TransducerElements // 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); // get image width and weight const unsigned int inputL = get_image_width( dSource ); const unsigned int inputS = get_image_height( dSource ) / (PAImage + 1); const unsigned int Slices = get_image_depth( dSource ); // 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 < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS; float tan_phi = tan(Angle / 360 * 2 * M_PI); - float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch; + float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch * outputL / TransducerElements; float part = part_multiplicator * s_i; if (part < 1) part = 1; short maxLine = min((l_i + part) + 1, (float)inputL); short minLine = max((l_i - part), 0.0f); short usedLines = (maxLine - minLine); float apod_mult = apodArraySize / (maxLine - minLine); short AddSample = 0; float output = 0; - float delayMultiplicator = pow(inputS / (RecordTime*SpeedOfSound) * Pitch, 2) / s_i / 2; + float delayMultiplicator = pow(inputS / (RecordTime*SpeedOfSound) * Pitch * 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 += apodArray[(short)((l_s - minLine)*apod_mult)] * read_imagef( dSource, defaultSampler, (int4)(l_s, AddSample, globalPosZ, 0 )).x; else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / usedLines; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl b/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl index 38509133f8..d8980e57fe 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl @@ -1,80 +1,81 @@ /*=================================================================== 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 ckDASSphe( __read_only image3d_t dSource, // input image __global float* dDest, // output buffer __global float* apodArray, unsigned short apodArraySize, unsigned short outputS, unsigned short outputL, float SpeedOfSound, float RecordTime, float Pitch, float Angle, - unsigned short PAImage // parameters + unsigned short PAImage, + unsigned short TransducerElements // 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); // get image width and weight const unsigned int inputL = get_image_width( dSource ); const unsigned int inputS = get_image_height( dSource ) / (PAImage + 1); const unsigned int Slices = get_image_depth( dSource ); // 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 < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS; float tan_phi = tan(Angle / 360 * 2 * M_PI); - float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch; + float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch * inputL / TransducerElements; float part = part_multiplicator * s_i; if (part < 1) part = 1; short maxLine = min((l_i + part) + 1, (float)inputL); short minLine = max((l_i - part), 0.0f); short usedLines = (maxLine - minLine); float apod_mult = apodArraySize / (maxLine - minLine); short AddSample = 0; float output = 0; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = sqrt( pow(s_i, 2) + - pow((inputS / (RecordTime*SpeedOfSound) * ((l_s - l_i)*Pitch*outputL) / inputL), 2) + pow((inputS / (RecordTime*SpeedOfSound) * ((l_s - l_i)*Pitch*TransducerElements) / inputL), 2) ); if (AddSample < inputS && AddSample >= 0) output += apodArray[(short)((l_s - minLine)*apod_mult)] * read_imagef( dSource, defaultSampler, (int4)(l_s, AddSample, globalPosZ, 0 )).x; else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / usedLines; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl index 45a7812795..a54dd151d5 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl @@ -1,125 +1,126 @@ /*=================================================================== 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 ckDMASQuad( __read_only image3d_t dSource, // input image __global float* dDest, // output buffer __global float* apodArray, unsigned short apodArraySize, unsigned short outputS, unsigned short outputL, float SpeedOfSound, float RecordTime, float Pitch, float Angle, - unsigned short PAImage // parameters + unsigned short PAImage, + unsigned short TransducerElements // 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); // get image width and weight const unsigned int inputL = get_image_width( dSource ); const unsigned int inputS = get_image_height( dSource ) / (PAImage + 1); const unsigned int Slices = get_image_depth( dSource ); // 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 < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS; float tan_phi = tan(Angle / 360 * 2 * M_PI); - float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch; + float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch * outputL / TransducerElements; float part = part_multiplicator * s_i; if (part < 1) part = 1; short maxLine = min((l_i + part) + 1, (float)inputL); short minLine = max((l_i - part), 0.0f); short usedLines = (maxLine - minLine); float apod_mult = apodArraySize / (maxLine - minLine); - float delayMultiplicator = pow((inputS / (RecordTime*SpeedOfSound) * Pitch), 2) / s_i / 2; + float delayMultiplicator = pow((inputS / (RecordTime*SpeedOfSound) * Pitch * TransducerElements / inputL), 2) / s_i / 2; float mult = 0; float output = 0; float AddSample1 = 0; float AddSample2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { - //AddSample1 = delayMultiplicator * pow((l_s1 - l_i), 2) + s_i; + AddSample1 = delayMultiplicator * pow((l_s1 - l_i), 2) + s_i; if (AddSample1 < inputS && AddSample1 >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { - //AddSample2 = delayMultiplicator * pow((l_s2 - l_i), 2) + s_i; + AddSample2 = delayMultiplicator * pow((l_s2 - l_i), 2) + s_i; if (AddSample2 < inputS && AddSample2 >= 0) { mult = read_imagef( dSource, defaultSampler, (int4)(l_s2, AddSample2, globalPosZ, 0 )).x * apodArray[(short)((l_s2 - minLine)*apod_mult)] * read_imagef( dSource, defaultSampler, (int4)(l_s1, AddSample1, globalPosZ, 0 )).x * apodArray[(short)((l_s1 - minLine)*apod_mult)]; - output += sqrt(mult*(mult>0)) * ((mult > 0) - (mult < 0)); + output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = 10 * output / (pow((float)usedLines, 2.0f) - (usedLines - 1)); } } /* //calculate the AddSamples beforehand to save some time for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i; } float mult = 0; float output = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { mult = read_imagef( dSource, defaultSampler, (int4)(l_s2, AddSample[l_s2], globalPosZ, 0 )).x * apodArray[(short)((l_s2 - minLine)*apod_mult)] * read_imagef( dSource, defaultSampler, (int4)(l_s1, AddSample[l_s1], globalPosZ, 0 )).x * apodArray[(short)((l_s1 - minLine)*apod_mult)]; output += sqrt(mult*(mult>0)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } */ \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl index aba1c30b68..9487abc201 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl @@ -1,97 +1,98 @@ /*=================================================================== 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 ckDMASSphe( __read_only image3d_t dSource, // input image __global float* dDest, // output buffer __global float* apodArray, unsigned short apodArraySize, unsigned short outputS, unsigned short outputL, float SpeedOfSound, float RecordTime, float Pitch, float Angle, - unsigned short PAImage // parameters + unsigned short PAImage, + unsigned short TransducerElements // 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); // get image width and weight const unsigned int inputL = get_image_width( dSource ); const unsigned int inputS = get_image_height( dSource ) / (PAImage + 1); const unsigned int Slices = get_image_depth( dSource ); // 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 < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS; float tan_phi = tan(Angle / 360 * 2 * M_PI); - float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch; + float part_multiplicator = tan_phi * RecordTime / inputS * SpeedOfSound / Pitch * outputL / TransducerElements; float part = part_multiplicator * s_i; if (part < 1) part = 1; short maxLine = min((l_i + part) + 1, (float)inputL); short minLine = max((l_i - part), 0.0f); short usedLines = (maxLine - minLine); float apod_mult = apodArraySize / (maxLine - minLine); float mult = 0; float output = 0; float AddSample1 = 0; float AddSample2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { AddSample1 = sqrt( pow(s_i, 2) + - pow((inputS / (RecordTime*SpeedOfSound) * ((l_s1 - l_i)*Pitch)), 2)); + pow((inputS / (RecordTime*SpeedOfSound) * ((l_s1 - l_i)*Pitch*TransducerElements)/inputL), 2)); if (AddSample1 < inputS && AddSample1 >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { AddSample2 = sqrt( pow(s_i, 2) + - pow((inputS / (RecordTime*SpeedOfSound) * ((l_s2 - l_i)*Pitch)), 2)); + pow((inputS / (RecordTime*SpeedOfSound) * ((l_s2 - l_i)*Pitch*TransducerElements)/inputL), 2)); if (AddSample2 < inputS && AddSample2 >= 0) { mult = read_imagef( dSource, defaultSampler, (int4)(l_s2, AddSample2, globalPosZ, 0 )).x * apodArray[(short)((l_s2 - minLine)*apod_mult)] * read_imagef( dSource, defaultSampler, (int4)(l_s1, AddSample1, globalPosZ, 0 )).x * apodArray[(short)((l_s1 - minLine)*apod_mult)]; - output += sqrt(mult*(mult>0)) * ((mult > 0) - (mult < 0)); + output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = 10 * output / (pow((float)usedLines, 2.0f) - (usedLines - 1)); - }*/ + } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/PixelSum.cl b/Modules/PhotoacousticsAlgorithms/Resources/PixelSum.cl new file mode 100644 index 0000000000..cbef559f47 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/PixelSum.cl @@ -0,0 +1,47 @@ +/*=================================================================== + +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 ckSum( + __read_only image3d_t dSource, // input image + __global float* dDest, // output buffer + __global float* input, + unsigned short num_images // 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); + + // get image width and weight + const unsigned int inputL = get_image_width( dSource ); + const unsigned int inputS = get_image_height( dSource ); + const unsigned int Slices = get_image_depth( dSource ); + + // 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 < inputL && globalPosY < inputS && globalPosZ < Slices ) + { + float output = 0; + for (unsigned short im = 0; im < num_images; ++im) + { + output += input[ im * Slices * inputL * inputS + globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX]; + } + dDest[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ] = output; + } +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index 81960ec7d9..c5a23a7339 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,14 +1,15 @@ set(CPP_FILES mitkPhotoacousticImage.cpp Algorithms/mitkPhotoacousticBeamformingFilter.cpp Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp ) set(RESOURCE_FILES DASQuadratic.cl DMASQuadratic.cl DASspherical.cl DMASspherical.cl + PixelSum.cl ) \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp index 8b34112494..8046d98378 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp @@ -1,265 +1,265 @@ /*=================================================================== 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 "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::Image< float, 3 > itkFloatImageType; typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); typedef itk::LinearInterpolateImageFunction T_Interpolator; itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); itkFloatImageType::SpacingType 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, int minSlice, int maxSlice) { float* inputData = new float[inputImage->GetDimension(0)*inputImage->GetDimension(1)*inputImage->GetDimension(2)]; float* outputData = new float[(inputImage->GetDimension(0) - left - right)*(inputImage->GetDimension(1) - above - below)*(maxSlice - minSlice + 1)]; 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, maxSlice - minSlice + 1 }; ImageReadAccessor acc(inputImage); // convert the data to float by default // as of now only those float, short, float are used at all... though it's easy to add other ones if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") { inputData = (float*)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]] = (float)inputPuffer[l*inputDim[1] + s + sl*inputDim[0]*inputDim[1]]; } } } } - else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") + else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (double)") { double* inputPuffer = (double*)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]] = (float)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 + minSlice)*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::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"); if (!config.partial) { config.bounds[0] = 0; config.bounds[1] = inputImage->GetDimension(2) - 1; } Image::Pointer croppedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0, config.bounds[0], config.bounds[1]); 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 BeamformingFilter::Pointer Beamformer = BeamformingFilter::New(); Beamformer->SetInput(resizedImage); Beamformer->Configure(config); Beamformer->SetProgressHandle(progressHandle); Beamformer->UpdateLargestPossibleRegion(); return Beamformer->GetOutput(); } \ No newline at end of file