diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp index 5aaf58f4ee..a817da318a 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDASFilter.cpp @@ -1,557 +1,558 @@ /*=================================================================== 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_Conf.Pitch = 0.0003; m_Conf.SpeedOfSound = 1540; m_Conf.SamplesPerLine = 2048; m_Conf.ReconstructionLines = 128; m_Conf.RecordTime = 0.00006; m_Conf.TransducerElements = 128; } 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()"); 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::BeamformingDASFilter::GenerateData() { + GenerateOutputInformation(); mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); double inputDim[2] = { input->GetDimension(0), input->GetDimension(1) / ((int)m_Conf.Photoacoustic + 1) }; double outputDim[2] = { output->GetDimension(0), output->GetDimension(1) }; const int apodArraySize = m_Conf.TransducerElements * 4; double* VonHannWindow = VonHannFunction(apodArraySize); if (!output->IsInitialized()) { return; } 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 double[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; m_InputDataPuffer = new double[input->GetDimension(0)*input->GetDimension(1)]; 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*(unsigned short)inputDim[1] + s] = (double)InputPuffer[l*(unsigned 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*(unsigned short)inputDim[1] + s] = (double)InputPuffer[l*(unsigned short)inputDim[1] + s]; } } m_InputData = m_InputDataPuffer; } else { MITK_INFO << "Could not determine pixel type"; return; } for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(unsigned short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(unsigned short)outputDim[0]]; for (unsigned short line = 0; line < outputDim[0]; ++line) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Linear) { threads[line] = std::thread(&BeamformingDASFilter::DASLinearLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { threads[line] = std::thread(&BeamformingDASFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { threads[line] = std::thread(&BeamformingDASFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); } } for (unsigned short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } output->SetSlice(m_OutputData, i); if (i % 50 == 0) MITK_INFO << "slice " << i + 1 << " of " << output->GetDimension(2) << " computed"; delete[] m_OutputData; delete[] m_InputDataPuffer; m_OutputData = nullptr; m_InputData = nullptr; } if (m_Conf.UseBP) { 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* VonHannWindow = new double[samples]; for (int n = 0; n < samples; ++n) { VonHannWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; } return VonHannWindow; } void mitk::BeamformingDASFilter::DASLinearLine(double* input, double* output, double inputDim[2], double outputDim[2], const unsigned short& line, double* apodisation, const unsigned short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; unsigned short AddSample = 0; unsigned short maxLine = 0; unsigned 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; + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double VH_mult = 1; double mult = 0; //linear delay l_i = line / outputL * inputL; l = (inputL / 2 - l_i) / inputL*m_Conf.Pitch*m_Conf.TransducerElements; for (unsigned short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (unsigned short)std::min((l_i + part) + 1, inputL); minLine = (unsigned short)std::max((l_i - part), 0.0); VH_mult = apodArraySize / (maxLine - minLine); x = m_Conf.RecordTime / inputS * s_i * m_Conf.SpeedOfSound; root = l / sqrt(pow(l, 2) + pow(m_Conf.RecordTime / inputS * s_i * m_Conf.SpeedOfSound, 2)); delayMultiplicator = root / (m_Conf.RecordTime*m_Conf.SpeedOfSound) *m_Conf.Pitch*m_Conf.TransducerElements / inputL; for (unsigned short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = delayMultiplicator * (l_s - l_i) + s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(unsigned short)outputL + line] += input[l_s + AddSample*(unsigned short)inputL] * apodisation[(unsigned short)((l_s - minLine)*VH_mult)]; } output[sample*(unsigned short)outputL + line] = output[sample*(unsigned short)outputL + line] / (maxLine - minLine); } } void mitk::BeamformingDASFilter::DASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const unsigned short& line, double* apodisation, const unsigned short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; unsigned short AddSample = 0; unsigned short maxLine = 0; unsigned 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; + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double VH_mult = 1; double mult = 0; //quadratic delay l_i = line / outputL * inputL; for (unsigned short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (unsigned short)std::min((l_i + part) + 1, inputL); minLine = (unsigned short)std::max((l_i - part), 0.0); VH_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 (unsigned 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*(unsigned short)outputL + line] += input[l_s + AddSample*(unsigned short)inputL] * apodisation[(unsigned short)((l_s - minLine)*VH_mult)]; } } output[sample*(unsigned short)outputL + line] = output[sample*(unsigned short)outputL + line] / (maxLine - minLine); } } void mitk::BeamformingDASFilter::DASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const unsigned short& line, double* apodisation, const unsigned short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; unsigned short AddSample = 0; unsigned short maxLine = 0; unsigned 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; + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double VH_mult = 1; double mult = 0; //exact delay l_i = (double)line / outputL * inputL; for (unsigned short sample = 0; sample < outputS; ++sample) { s_i = (double)sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (unsigned short)std::min((l_i + part) + 1, inputL); minLine = (unsigned short)std::max((l_i - part), 0.0); VH_mult = apodArraySize / (maxLine - minLine); for (unsigned 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*(unsigned short)outputL + line] += input[l_s + AddSample*(unsigned short)inputL] * apodisation[(unsigned short)((l_s - minLine)*VH_mult)]; } } output[sample*(unsigned short)outputL + line] = output[sample*(unsigned short)outputL + line] / (maxLine - minLine); } } mitk::Image::Pointer mitk::BeamformingDASFilter::BandpassFilter(mitk::Image::Pointer data) { 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;*/ //debugging 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(toReal->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) { // tukey window double alpha = m_Conf.BPFalloff; 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; } 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; } } 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.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp index bac4d1d585..5d471af217 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingDMASFilter.cpp @@ -1,655 +1,655 @@ /*=================================================================== 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 "mitkProperties.h" #include "mitkImageReadAccessor.h" #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) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_Conf.Pitch = 0.0003; m_Conf.SpeedOfSound = 1540; m_Conf.SamplesPerLine = 2048; m_Conf.ReconstructionLines = 128; m_Conf.RecordTime = 0.00006; m_Conf.TransducerElements = 128; } mitk::BeamformingDMASFilter::~BeamformingDMASFilter() { } void mitk::BeamformingDMASFilter::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() { 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::BeamformingDMASFilter::GenerateData() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); double inputDim[2] = { input->GetDimension(0), input->GetDimension(1) / ((int)m_Conf.Photoacoustic + 1) }; double outputDim[2] = { output->GetDimension(0), output->GetDimension(1) }; const int apodArraySize = m_Conf.TransducerElements * 4; double* VonHannWindow = VonHannFunction(apodArraySize); if (!output->IsInitialized()) { return; } 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 double[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; m_InputDataPuffer = new double[input->GetDimension(0)*input->GetDimension(1)]; 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*(unsigned short)inputDim[1] + s] = (double)InputPuffer[l*(unsigned 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*(unsigned short)inputDim[1] + s] = (double)InputPuffer[l*(unsigned short)inputDim[1] + s]; } } m_InputData = m_InputDataPuffer; } else { MITK_INFO << "Could not determine pixel type"; return; } for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(unsigned short)outputDim[1] + s] = 0; } } const unsigned short max_threads = 12; std::thread *threads = new std::thread[(unsigned short)outputDim[0]]; for (unsigned short line = 0; line < outputDim[0]; ++line) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Linear) { threads[line] = std::thread(&BeamformingDMASFilter::DMASLinearLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { threads[line] = std::thread(&BeamformingDMASFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { threads[line] = std::thread(&BeamformingDMASFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); } } for (unsigned short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } /* bool *threadfinished = new bool[max_threads]; for (unsigned short r = 0; r < max_threads; ++r) { threadfinished[r] = false; } unsigned short line = 0; while(line < outputDim[0]) { //DMASSphericalLine(m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize) //threads[line] = std::thread(&BeamformingDMASFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize); for (unsigned short n = 0; n < max_threads; ++n) { if (threadfinished[n]) { threads[n].join(); threadfinished[n] = false; //MITK_INFO << "thread " << n << " joined"; } if (!threads[n].joinable()) { threads[n] = std::thread(&BeamformingDMASFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, VonHannWindow, apodArraySize, &threadfinished[n]); ++line; //MITK_INFO << "thread " << n << " created for line " << line - 1; break; } } } for (unsigned short n = 0; n < max_threads; ++n) { if (threads[n].joinable()) { threads[n].join(); threadfinished[n] = false; //MITK_INFO << "thread " << n << " joined"; } } delete[] threadfinished;*/ //threadpool... seems slower output->SetSlice(m_OutputData, i); if (i % 50 == 0) MITK_INFO << "slice " << i + 1 << " of " << output->GetDimension(2) << " computed"; delete[] m_OutputData; delete[] m_InputDataPuffer; delete[] threads; m_OutputData = nullptr; m_InputData = nullptr; } if (m_Conf.UseBP) { 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[] VonHannWindow; 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; } void mitk::BeamformingDMASFilter::Configure(beamformingSettings settings) { m_Conf = settings; } mitk::Image::Pointer mitk::BeamformingDMASFilter::BandpassFilter(mitk::Image::Pointer data) { 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;*/ //debugging 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(toReal->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) { // tukey window double alpha = m_Conf.BPFalloff; 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; } 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; } } 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* VonHannWindow = new double[samples]; for (int n = 0; n < samples; ++n) { VonHannWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; } return VonHannWindow; } void mitk::BeamformingDMASFilter::DMASLinearLine(double* input, double* output, double inputDim[2], double outputDim[2], const unsigned short& line, double* apodisation, const unsigned short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; unsigned short AddSample1 = 0; unsigned short AddSample2 = 0; unsigned short maxLine = 0; unsigned 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; + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double VH_mult = 1; double mult = 0; //linear delay l_i = line / outputL * inputL; l = (inputL / 2 - l_i) / inputL*m_Conf.Pitch*m_Conf.TransducerElements; for (unsigned short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (unsigned short)std::min((l_i + part) + 1, inputL); minLine = (unsigned short)std::max((l_i - part), 0.0); VH_mult = apodArraySize / (maxLine - minLine); x = m_Conf.RecordTime / inputS * s_i * m_Conf.SpeedOfSound; root = l / sqrt(pow(l, 2) + pow(m_Conf.RecordTime / inputS * s_i * m_Conf.SpeedOfSound, 2)); delayMultiplicator = root / (m_Conf.RecordTime*m_Conf.SpeedOfSound) *m_Conf.Pitch*m_Conf.TransducerElements / inputL; //calculate the AddSamples beforehand to save some time unsigned short* AddSample = new unsigned short[maxLine - minLine]; for (unsigned short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (unsigned short)(delayMultiplicator * (minLine + l_s - l_i) + s_i); } for (unsigned short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (unsigned 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] * (unsigned short)inputL] * input[l_s1 + AddSample[l_s1 - minLine] * (unsigned short)inputL] * apodisation[(unsigned short)((l_s1 - minLine)*VH_mult)] * apodisation[(unsigned short)((l_s2 - minLine)*VH_mult)]; output[sample*(unsigned short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } } output[sample*(unsigned short)outputL + line] = 10 * output[sample*(unsigned short)outputL + line] / (pow(maxLine - minLine, 2) - (maxLine - minLine)); delete[] AddSample; } } void mitk::BeamformingDMASFilter::DMASQuadraticLine(double* input, double* output, double inputDim[2], double outputDim[2], const unsigned short& line, double* apodisation, const unsigned short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; unsigned short AddSample1 = 0; unsigned short AddSample2 = 0; unsigned short maxLine = 0; unsigned 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; + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double VH_mult = 1; double mult = 0; //quadratic delay l_i = line / outputL * inputL; for (unsigned short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (unsigned short)std::min((l_i + part) + 1, inputL); minLine = (unsigned short)std::max((l_i - part), 0.0); VH_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 unsigned short* AddSample = new unsigned short[maxLine - minLine]; for (unsigned short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (unsigned short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i); } for (unsigned short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (unsigned 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] * (unsigned short)inputL] * apodisation[(unsigned short)((l_s2 - minLine)*VH_mult)] * input[l_s1 + AddSample[l_s1 - minLine] * (unsigned short)inputL] * apodisation[(unsigned short)((l_s1 - minLine)*VH_mult)]; output[sample*(unsigned short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } } output[sample*(unsigned short)outputL + line] = 10 * output[sample*(unsigned short)outputL + line] / (pow(maxLine - minLine, 2) - (maxLine - minLine - 1)); delete[] AddSample; } } void mitk::BeamformingDMASFilter::DMASSphericalLine(double* input, double* output, double inputDim[2], double outputDim[2], const unsigned short& line, double* apodisation, const unsigned short& apodArraySize) { double& inputS = inputDim[1]; double& inputL = inputDim[0]; double& outputS = outputDim[1]; double& outputL = outputDim[0]; unsigned short AddSample1 = 0; unsigned short AddSample2 = 0; unsigned short maxLine = 0; unsigned 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; + double part_multiplicator = tan_phi * m_Conf.RecordTime / inputS * m_Conf.SpeedOfSound / m_Conf.Pitch * m_Conf.ReconstructionLines / m_Conf.TransducerElements; double VH_mult = 1; double mult = 0; //exact delay l_i = line / outputL * inputL; for (unsigned short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (unsigned short)std::min((l_i + part) + 1, inputL); minLine = (unsigned short)std::max((l_i - part), 0.0); VH_mult = apodArraySize / (maxLine - minLine); //calculate the AddSamples beforehand to save some time unsigned short* AddSample = new unsigned short[maxLine - minLine]; for (unsigned short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (unsigned 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 (unsigned short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (unsigned 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] * (unsigned short)inputL] * apodisation[(unsigned short)((l_s2 - minLine)*VH_mult)] * input[l_s1 + AddSample[l_s1 - minLine] * (unsigned short)inputL] * apodisation[(unsigned short)((l_s1 - minLine)*VH_mult)]; output[sample*(unsigned short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } } output[sample*(unsigned short)outputL + line] = 10 * output[sample*(unsigned short)outputL + line] / (pow(maxLine - minLine, 2) - (maxLine - minLine - 1)); delete[] AddSample; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt index 4fc45dbea1..4a47e80aba 100644 --- a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt +++ b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt @@ -1,8 +1,8 @@ MITK_CREATE_MODULE( SUBPROJECTS - DEPENDS MitkCore + DEPENDS MitkCore MitkAlgorithmsExt #AUTOLOAD_WITH MitkCore INCLUDE_DIRS PUBLIC Algorithms/ITKUltrasound INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} PACKAGE_DEPENDS ITK|ITKFFT+ITKImageCompose+ITKImageIntensity ) diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp index 83685afaa7..6e790ceabc 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp @@ -1,202 +1,224 @@ /*=================================================================== 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 "Algorithms/ITKUltrasound/itkBModeImageFilter.h" #include "Algorithms/itkPhotoacousticBModeImageFilter.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" #include "Algorithms/mitkPhotoacousticBeamformingDASFilter.h" #include "Algorithms/mitkPhotoacousticBeamformingDMASFilter.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::ApplyBeamformingDAS(mitk::Image::Pointer inputImage, BeamformingDASFilter::beamformingSettings config) +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left) { - BeamformingDASFilter::Pointer Beamformer = BeamformingDASFilter::New(); + mitk::AutoCropImageFilter::Pointer cropFilter = mitk::AutoCropImageFilter::New(); + + itk::ImageRegion<3> cropRegion; + itk::Index<3> index = { left, above, 0 }; + itk::Size<3> size = { inputImage->GetDimension(0) - left - right, inputImage->GetDimension(1) - above - below, inputImage->GetDimension(2) }; + cropRegion.SetIndex(index); + cropRegion.SetSize(size); + + cropFilter->SetInput(inputImage); + cropFilter->SetCroppingRegion(cropRegion); + cropFilter->Update(); - Image::Pointer resizedImage = inputImage; + return cropFilter->GetOutput(); +} + +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamformingDAS(mitk::Image::Pointer inputImage, BeamformingDASFilter::beamformingSettings config, int cutoff) +{ + config.RecordTime = config.RecordTime - cutoff / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping + Image::Pointer croppedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0); + Image::Pointer resizedImage = croppedImage; if (inputImage->GetDimension(0) != config.ReconstructionLines) { auto begin = std::chrono::high_resolution_clock::now(); unsigned int dim[3] = { config.ReconstructionLines, inputImage->GetDimension(1), inputImage->GetDimension(2) }; resizedImage = ApplyResampling(inputImage, 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; } + BeamformingDASFilter::Pointer Beamformer = BeamformingDASFilter::New(); + Beamformer->SetInput(resizedImage); Beamformer->Configure(config); Beamformer->UpdateLargestPossibleRegion(); return Beamformer->GetOutput(); } -mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamformingDMAS(mitk::Image::Pointer inputImage, BeamformingDMASFilter::beamformingSettings config) +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamformingDMAS(mitk::Image::Pointer inputImage, BeamformingDMASFilter::beamformingSettings config, int cutoff) { - BeamformingDMASFilter::Pointer Beamformer = BeamformingDMASFilter::New(); - - Image::Pointer resizedImage = inputImage; + config.RecordTime = config.RecordTime - cutoff / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping + Image::Pointer croppedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0); + Image::Pointer resizedImage = croppedImage; if(inputImage->GetDimension(0) != config.ReconstructionLines) { auto begin = std::chrono::high_resolution_clock::now(); unsigned int dim[3] = { config.ReconstructionLines, inputImage->GetDimension(1), inputImage->GetDimension(2) }; resizedImage = ApplyResampling(inputImage, 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; } + BeamformingDMASFilter::Pointer Beamformer = BeamformingDMASFilter::New(); + Beamformer->SetInput(resizedImage); Beamformer->Configure(config); Beamformer->UpdateLargestPossibleRegion(); return Beamformer->GetOutput(); } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h index b916962f77..188423fa84 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h @@ -1,49 +1,50 @@ /*=================================================================== 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 "Algorithms/mitkPhotoacousticBeamformingDASFilter.h" #include "Algorithms/mitkPhotoacousticBeamformingDMASFilter.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); - mitk::Image::Pointer ApplyBeamformingDMAS(mitk::Image::Pointer inputImage, BeamformingDMASFilter::beamformingSettings config); + mitk::Image::Pointer ApplyBeamformingDAS(mitk::Image::Pointer inputImage, BeamformingDASFilter::beamformingSettings config, int cutoff); + mitk::Image::Pointer ApplyBeamformingDMAS(mitk::Image::Pointer inputImage, BeamformingDMASFilter::beamformingSettings config, int cutoff); + 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 930cb5058e..a70a633ac2 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,379 +1,379 @@ /*=================================================================== 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 //mitk image #include #include "mitkPhotoacousticImage.h" #include "Algorithms\mitkPhotoacousticBeamformingDASFilter.h" #include "Algorithms\mitkPhotoacousticBeamformingDMASFilter.h" const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false) { } 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(ApplyBModeFilter())); 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(ApplyBeamforming())); 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(double)), this, SLOT(UpdateFrequency())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked), this, SLOT(UpdateFrequency())); m_Controls.DoResampling->setChecked(false); m_Controls.ResamplingValue->setEnabled(false); UseImageSpacing(); } void PAImageProcessing::UpdateFrequency() { DMASconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] DMASconfig.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 (image) UpdateRecordTime(image); } std::stringstream frequency; frequency << 1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000; frequency << "Hz"; m_Controls.BPInfoDisplay->setText(frequency.str().c_str()); } 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::ApplyBModeFilter() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information( NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { std::stringstream message; std::string name; message << "Performing image processing for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; } message << "."; MITK_INFO << message.str(); auto newNode = mitk::DataNode::New(); mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); newNode->SetData(filterbank->ApplyBmodeFilter(image, m_UseLogfilter, m_ResampleSpacing)); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data),true,true); newNode->SetLevelWindow(levelWindow); // name the new Data node std::stringstream newNodeName; if (node->GetName(name)) { // a property called "name" was found for this DataNode newNodeName << name << ", "; } newNodeName << "B-Mode"; newNode->SetName(newNodeName.str()); // add new node to data storage this->GetDataStorage()->Add(newNode); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } } void PAImageProcessing::ApplyBeamforming() { 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 << "'"; } message << "."; MITK_INFO << message.str(); mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); auto newNode = mitk::DataNode::New(); if(m_CurrentBeamformingAlgorithm == BeamformingAlgorithms::DAS) - newNode->SetData(filterbank->ApplyBeamformingDAS(image, DASconfig)); + newNode->SetData(filterbank->ApplyBeamformingDAS(image, DASconfig, m_Controls.Cutoff->value())); else if(m_CurrentBeamformingAlgorithm == BeamformingAlgorithms::DMAS) - newNode->SetData(filterbank->ApplyBeamformingDMAS(image, DMASconfig)); + newNode->SetData(filterbank->ApplyBeamformingDMAS(image, DMASconfig, m_Controls.Cutoff->value())); // name the new Data node std::stringstream newNodeName; if (node->GetName(name)) { // a property called "name" was found for this DataNode newNodeName << name << ", "; } if (m_CurrentBeamformingAlgorithm == BeamformingAlgorithms::DAS) newNodeName << "DAS beamformed, "; else if (m_CurrentBeamformingAlgorithm == BeamformingAlgorithms::DMAS) newNodeName << "DMAS beamformed, "; if (DASconfig.DelayCalculationMethod == mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::Linear) newNodeName << "linear delay"; if (DASconfig.DelayCalculationMethod == mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::QuadApprox) newNodeName << "quadratic delay"; if (DASconfig.DelayCalculationMethod == mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::Spherical) newNodeName << "spherical delay"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } } void PAImageProcessing::UpdateBFSettings(mitk::Image::Pointer image) { if ("Delay and Sum" == m_Controls.BFAlgorithm->currentText()) m_CurrentBeamformingAlgorithm = BeamformingAlgorithms::DAS; else if ("Delay, Multiply and Sum" == m_Controls.BFAlgorithm->currentText()) m_CurrentBeamformingAlgorithm = BeamformingAlgorithms::DMAS; if ("Linear Waves" == m_Controls.DelayCalculation->currentText()) { DASconfig.DelayCalculationMethod = mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::Linear; DMASconfig.DelayCalculationMethod = mitk::BeamformingDMASFilter::beamformingSettings::DelayCalc::Linear; } else if ("Quadratic Approximation to Spherical Waves" == m_Controls.DelayCalculation->currentText()) { DASconfig.DelayCalculationMethod = mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::QuadApprox; DMASconfig.DelayCalculationMethod = mitk::BeamformingDMASFilter::beamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Waves" == m_Controls.DelayCalculation->currentText()) { DASconfig.DelayCalculationMethod = mitk::BeamformingDASFilter::beamformingSettings::DelayCalc::Spherical; DMASconfig.DelayCalculationMethod = mitk::BeamformingDMASFilter::beamformingSettings::DelayCalc::Spherical; } 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(); DASconfig.BPHighPass = 1000000 * m_Controls.BPhigh->value(); DASconfig.BPLowPass = 1000000 * (1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000 / 1000000 - m_Controls.BPlow->value()); DASconfig.BPFalloff = 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(); DMASconfig.BPHighPass = 1000000 * m_Controls.BPhigh->value(); DMASconfig.BPLowPass = 1000000 * (1 / (DMASconfig.RecordTime / DMASconfig.SamplesPerLine) * DMASconfig.SamplesPerLine / 2 / 2 / 1000 /1000000 - m_Controls.BPlow->value()); DMASconfig.BPFalloff = m_Controls.BPFalloff->value(); DMASconfig.UseBP = m_Controls.UseBP->isChecked(); UpdateRecordTime(image); } 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"; } else { DASconfig.RecordTime = m_Controls.ScanDepth->value() / 1000 / DASconfig.SpeedOfSound; // [s] DMASconfig.RecordTime = m_Controls.ScanDepth->value() / 1000 / DMASconfig.SpeedOfSound; // [s] } if ("Ultrasound Image" == m_Controls.ImageType->currentText()) { if (m_Controls.UseImageSpacing->isChecked()) { - //DASconfig.RecordTime = DASconfig.RecordTime / 2; // [s] - //DMASconfig.RecordTime = DMASconfig.RecordTime / 2; // [s] + DASconfig.RecordTime = DASconfig.RecordTime / 2; // [s] + DMASconfig.RecordTime = DMASconfig.RecordTime / 2; // [s] } DASconfig.Photoacoustic = false; DMASconfig.Photoacoustic = false; } else { DASconfig.Photoacoustic = true; DMASconfig.Photoacoustic = true; } } void PAImageProcessing::UseImageSpacing() { if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ScanDepth->setDisabled(true); } else { m_Controls.ScanDepth->setEnabled(true); } } \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui index c70f85f4f7..0ea1b864b9 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui @@ -1,493 +1,513 @@ PAImageProcessingControls 0 0 615 678 0 0 QmitkTemplate <html><head/><body><p><span style=" font-size:10pt; font-weight:600;">B-mode Filter</span></p></body></html> <html><head/><body><p><span style=" font-weight:600;">Filter Settings</span></p></body></html> Do Resampling 0.010000000000000 1.000000000000000 0.010000000000000 0.150000000000000 Resample Spacing [mm] Add Logfilter QLabel { color: rgb(255, 0, 0) } <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600;">Please select an image!</span></p></body></html> Do image processing Apply B-mode Filter Qt::Horizontal <html><head/><body><p><span style=" font-size:10pt; font-weight:600;">Beamforming</span></p></body></html> <html><head/><body><p><span style=" font-weight:600;">Settings</span></p></body></html> - + Speed of Sound [m/s] - + Scan Depth [mm] - + 1 0.100000000000000 2.000000000000000 0.100000000000000 0.300000000000000 - + 1 200.000000000000000 3000.000000000000000 0.100000000000000 1540.000000000000000 - + Beamforming Method - + Transducer Element Count - + 1 0.100000000000000 50.000000000000000 - + Delay and Sum Delay, Multiply and Sum - + Transducer Pitch [mm] - + 1 1024 128 128 - + 256 16384 256 2048 - + 128 2048 128 128 - + Samples - + Reconstruction Lines - + Linear Waves Quadratic Approximation to Spherical Waves Spherical Waves - + Delay Calculation Method - + 1 1.000000000000000 10.000000000000000 - + Compute ScanDepth from Image Data true - + Element Angle [°] - + Photoacoustic Image Ultrasound Image - + Image Type - + 2 0.010000000000000 1.000000000000000 0.100000000000000 0.200000000000000 - + 2 20.000000000000000 - + 2 10.000000000000000 200.000000000000000 - 50.000000000000000 + 20.000000000000000 - + High Pass BP Frequency [MHz] - + BP Falloff (0: rectangular, 1: smooth) - + Low Pass BP Frequency [MHz] - + Frequency Range: - + - - + Use Bandpass After Beamforming true + + + + 500 + + + 10 + + + 50 + + + + + + + Cutoff Upper Voxels + + + Apply Beamforming Qt::Vertical 20 40