diff --git a/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h index 7b1332a8c7..e201d468a5 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h @@ -1,63 +1,53 @@ /*=================================================================== 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 MITKVESSELDRAWER_H #define MITKVESSELDRAWER_H #include "mitkVector.h" #include "mitkPAVesselMeanderStrategy.h" #include "mitkPAInSilicoTissueVolume.h" #include "mitkPAVector.h" #include "mitkPAVesselProperties.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT VesselDrawer : public itk::LightObject { public: mitkClassMacroItkParent(VesselDrawer, itk::LightObject); itkFactorylessNewMacro(Self); void DrawVesselInVolume( VesselProperties::Pointer properties, InSilicoTissueVolume::Pointer volume); protected: VesselDrawer(); virtual ~VesselDrawer(); private: }; - /** - * @brief Equal A function comparing two vessel drawer instances for beeing equal - * - * @param rightHandSide A vesseldrawer to be compared - * @param leftHandSide A vesseldrawer to be compared - * @param eps tolarence for comparison. You can use mitk::eps in most cases. - * @param verbose flag indicating if the user wants detailed console output or not. - * @return true, if all subsequent comparisons are true, false otherwise - */ - MITKPHOTOACOUSTICSLIB_EXPORT bool Equal(const VesselDrawer::Pointer leftHandSide, const VesselDrawer::Pointer rightHandSide, double eps, bool verbose); - } + } } #endif // MITKVESSELDRAWER_H diff --git a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp index 9c9e38f29b..eb8e036a2b 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp @@ -1,365 +1,360 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include mitk::pa::InSilicoTissueVolume::InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters) { { unsigned int xDim = parameters->GetXDim(); unsigned int yDim = parameters->GetYDim(); unsigned int zDim = parameters->GetZDim(); m_TDim = 4; unsigned int size = xDim * yDim * zDim; double* absorptionArray = new double[size]; double* scatteringArray = new double[size]; double* anisotropyArray = new double[size]; double* segmentationArray = new double[size]; for (unsigned int index = 0; index < size; index++) { absorptionArray[index] = parameters->GetBackgroundAbsorption(); scatteringArray[index] = parameters->GetBackgroundScattering(); anisotropyArray[index] = parameters->GetBackgroundAnisotropy(); segmentationArray[index] = SegmentationType::BACKGROUND; } m_AbsorptionVolume = Volume::New(absorptionArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); m_ScatteringVolume = Volume::New(scatteringArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); m_AnisotropyVolume = Volume::New(anisotropyArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); m_SegmentationVolume = Volume::New(segmentationArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); } m_TissueParameters = parameters; m_PropertyList = mitk::PropertyList::New(); UpdatePropertyList(); } void mitk::pa::InSilicoTissueVolume::UpdatePropertyList() { //Set properties AddIntProperty("mcflag", m_TissueParameters->GetMCflag()); AddIntProperty("launchflag", m_TissueParameters->GetMCLaunchflag()); AddIntProperty("boundaryflag", m_TissueParameters->GetMCBoundaryflag()); AddDoubleProperty("launchPointX", m_TissueParameters->GetMCLaunchPointX()); AddDoubleProperty("launchPointY", m_TissueParameters->GetMCLaunchPointY()); AddDoubleProperty("launchPointZ", m_TissueParameters->GetMCLaunchPointZ()); AddDoubleProperty("focusPointX", m_TissueParameters->GetMCFocusPointX()); AddDoubleProperty("focusPointY", m_TissueParameters->GetMCFocusPointY()); AddDoubleProperty("focusPointZ", m_TissueParameters->GetMCFocusPointZ()); AddDoubleProperty("trajectoryVectorX", m_TissueParameters->GetMCTrajectoryVectorX()); AddDoubleProperty("trajectoryVectorY", m_TissueParameters->GetMCTrajectoryVectorY()); AddDoubleProperty("trajectoryVectorZ", m_TissueParameters->GetMCTrajectoryVectorZ()); AddDoubleProperty("radius", m_TissueParameters->GetMCRadius()); AddDoubleProperty("waist", m_TissueParameters->GetMCWaist()); AddDoubleProperty("partialVolume", m_TissueParameters->GetDoPartialVolume()); AddDoubleProperty("standardTissueAbsorption", m_TissueParameters->GetBackgroundAbsorption()); AddDoubleProperty("standardTissueScattering", m_TissueParameters->GetBackgroundScattering()); AddDoubleProperty("standardTissueAnisotropy", m_TissueParameters->GetBackgroundAnisotropy()); AddDoubleProperty("airThickness", m_TissueParameters->GetAirThicknessInMillimeters()); AddDoubleProperty("skinThickness", m_TissueParameters->GetSkinThicknessInMillimeters()); } mitk::pa::InSilicoTissueVolume::InSilicoTissueVolume( Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList) { m_AbsorptionVolume = absorptionVolume; m_ScatteringVolume = scatteringVolume; m_AnisotropyVolume = anisotropyVolume; m_SegmentationVolume = segmentationVolume; m_TissueParameters = tissueParameters; m_PropertyList = propertyList; if (m_SegmentationVolume.IsNotNull()) m_TDim = 4; else m_TDim = 3; } double mitk::pa::InSilicoTissueVolume::GetSpacing() { return m_AbsorptionVolume->GetSpacing(); } void mitk::pa::InSilicoTissueVolume::SetSpacing(double spacing) { m_AbsorptionVolume->SetSpacing(spacing); m_ScatteringVolume->SetSpacing(spacing); m_AnisotropyVolume->SetSpacing(spacing); m_SegmentationVolume->SetSpacing(spacing); } void mitk::pa::InSilicoTissueVolume::AddDoubleProperty(std::string label, double value) { m_PropertyList->SetDoubleProperty(label.c_str(), value); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New(label)); } void mitk::pa::InSilicoTissueVolume::AddIntProperty(std::string label, int value) { m_PropertyList->SetIntProperty(label.c_str(), value); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New(label)); } mitk::Image::Pointer mitk::pa::InSilicoTissueVolume::ConvertToMitkImage() { mitk::Image::Pointer resultImage = mitk::Image::New(); mitk::PixelType TPixel = mitk::MakeScalarPixelType(); unsigned int* dimensionsOfImage = new unsigned int[4]; // Copy dimensions dimensionsOfImage[0] = m_TissueParameters->GetYDim(); dimensionsOfImage[1] = m_TissueParameters->GetXDim(); dimensionsOfImage[2] = m_TissueParameters->GetZDim(); dimensionsOfImage[3] = m_TDim; - MITK_INFO << "Initializing image..."; resultImage->Initialize(TPixel, 4, dimensionsOfImage, 1); - MITK_INFO << "Initializing image...[Done]"; mitk::Vector3D spacing; spacing.Fill(m_TissueParameters->GetVoxelSpacingInCentimeters()); resultImage->SetSpacing(spacing); - MITK_INFO << "Set Import Volumes..."; - //Copy memory, deal with cleaning up memory yourself! resultImage->SetImportVolume(m_AbsorptionVolume->GetData(), 0, 0, mitk::Image::CopyMemory); resultImage->SetImportVolume(m_ScatteringVolume->GetData(), 1, 0, mitk::Image::CopyMemory); resultImage->SetImportVolume(m_AnisotropyVolume->GetData(), 2, 0, mitk::Image::CopyMemory); resultImage->SetImportVolume(m_SegmentationVolume->GetData(), 3, 0, mitk::Image::CopyMemory); - MITK_INFO << "Set Import Volumes...[Done]"; resultImage->SetPropertyList(m_PropertyList); return resultImage; } mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueVolume::New( Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList) { InSilicoTissueVolume::Pointer smartPtr = new InSilicoTissueVolume( absorptionVolume, scatteringVolume, anisotropyVolume, segmentationVolume, tissueParameters, propertyList); smartPtr->UnRegister(); return smartPtr; } mitk::pa::InSilicoTissueVolume::~InSilicoTissueVolume() { m_AbsorptionVolume = nullptr; m_ScatteringVolume = nullptr; m_AnisotropyVolume = nullptr; m_SegmentationVolume = nullptr; m_TissueParameters = nullptr; m_PropertyList = nullptr; } void mitk::pa::InSilicoTissueVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy) { if (IsInsideVolume(x, y, z)) { m_AbsorptionVolume->SetData(absorption, x, y, z); m_ScatteringVolume->SetData(scattering, x, y, z); m_AnisotropyVolume->SetData(anisotropy, x, y, z); } } void mitk::pa::InSilicoTissueVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType) { if (IsInsideVolume(x, y, z)) { m_AbsorptionVolume->SetData(absorption, x, y, z); m_ScatteringVolume->SetData(scattering, x, y, z); m_AnisotropyVolume->SetData(anisotropy, x, y, z); m_SegmentationVolume->SetData(segmentType, x, y, z); } } bool mitk::pa::InSilicoTissueVolume::IsInsideVolume(int x, int y, int z) { return x >= 0 && x < m_TissueParameters->GetXDim() && y >= 0 && y < m_TissueParameters->GetYDim() && z >= 0 && z < m_TissueParameters->GetZDim(); } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetAbsorptionVolume() { return m_AbsorptionVolume; } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetSegmentationVolume() { return m_SegmentationVolume; } void mitk::pa::InSilicoTissueVolume::FinalizeVolume() { AddSkinAndAirLayers(); // If specified, randomize all tissue parameters if (m_TissueParameters->GetRandomizePhysicalProperties()) RandomizeTissueCoefficients(m_TissueParameters->GetUseRngSeed(), m_TissueParameters->GetRngSeed(), m_TissueParameters->GetRandomizePhysicalPropertiesPercentage()); } void mitk::pa::InSilicoTissueVolume::AddSkinAndAirLayers() { //Calculate the index location according to thickness in cm double airvoxel = (m_TissueParameters->GetAirThicknessInMillimeters() / m_TissueParameters->GetVoxelSpacingInCentimeters()) / 10; double skinvoxel = airvoxel + (m_TissueParameters->GetSkinThicknessInMillimeters() / m_TissueParameters->GetVoxelSpacingInCentimeters()) / 10; for (int y = 0; y < m_TissueParameters->GetYDim(); y++) { for (int x = 0; x < m_TissueParameters->GetXDim(); x++) { // Add air from index 0 to airvoxel if (m_TissueParameters->GetAirThicknessInMillimeters() > mitk::eps) { FillZLayer(x, y, 0, airvoxel, m_TissueParameters->GetAirAbsorption(), m_TissueParameters->GetAirScattering(), m_TissueParameters->GetAirAnisotropy(), SegmentationType::AIR); } //Add skin from index airvoxel to skinvoxel if (m_TissueParameters->GetSkinThicknessInMillimeters() > mitk::eps) { FillZLayer(x, y, airvoxel, skinvoxel, m_TissueParameters->GetSkinAbsorption(), m_TissueParameters->GetSkinScattering(), m_TissueParameters->GetSkinAnisotropy(), SegmentationType::SKIN); } } } } void mitk::pa::InSilicoTissueVolume::FillZLayer(int x, int y, double startIdx, double endIdx, double absorption, double scattering, double anisotropy, SegmentationType segmentationType) { for (int z = startIdx; z < endIdx; z++) { if (IsInsideVolume(x, y, z)) { if (endIdx - z < 1) { //Simulate partial volume effects m_AbsorptionVolume->SetData((1 - (endIdx - z)) * m_AbsorptionVolume->GetData(x, y, z) + (endIdx - z) * absorption, x, y, z); m_ScatteringVolume->SetData((1 - (endIdx - z)) * m_ScatteringVolume->GetData(x, y, z) + (endIdx - z) * scattering, x, y, z); m_AnisotropyVolume->SetData((1 - (endIdx - z)) * m_AnisotropyVolume->GetData(x, y, z) + (endIdx - z) * anisotropy, x, y, z); if (endIdx - z > 0.5) { //Only put the segmentation label if more than half of the partial volume is the wanted tissue type m_SegmentationVolume->SetData(segmentationType, x, y, z); } } else { m_AbsorptionVolume->SetData(absorption, x, y, z); m_ScatteringVolume->SetData(scattering, x, y, z); m_AnisotropyVolume->SetData(anisotropy, x, y, z); m_SegmentationVolume->SetData(segmentationType, x, y, z); } } } } void mitk::pa::InSilicoTissueVolume::RandomizeTissueCoefficients(long rngSeed, bool useRngSeed, double percentage) { std::mt19937 rng; std::random_device randomDevice; if (useRngSeed) { rng.seed(rngSeed); } else { if (randomDevice.entropy() > 0.1) { rng.seed(randomDevice()); } else { rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); } } std::normal_distribution<> percentageDistribution(1, percentage / 100); for (int y = 0; y < m_TissueParameters->GetYDim(); y++) { for (int x = 0; x < m_TissueParameters->GetXDim(); x++) { for (int z = 0; z < m_TissueParameters->GetZDim(); z++) { m_AbsorptionVolume->SetData(m_AbsorptionVolume->GetData(x, y, z) * percentageDistribution(rng), x, y, z); m_ScatteringVolume->SetData(m_ScatteringVolume->GetData(x, y, z) * percentageDistribution(rng), x, y, z); } } } } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetScatteringVolume() { return m_ScatteringVolume; } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetAnisotropyVolume() { return m_AnisotropyVolume; } void mitk::pa::InSilicoTissueVolume::SetAbsorptionVolume(Volume::Pointer volume) { m_AbsorptionVolume = volume; } void mitk::pa::InSilicoTissueVolume::SetScatteringVolume(Volume::Pointer volume) { m_ScatteringVolume = volume; } void mitk::pa::InSilicoTissueVolume::SetAnisotropyVolume(Volume::Pointer volume) { m_AnisotropyVolume = volume; } void mitk::pa::InSilicoTissueVolume::SetSegmentationVolume(Volume::Pointer volume) { m_SegmentationVolume = volume; } diff --git a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp index fd65f072c4..992a9ce378 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp @@ -1,150 +1,153 @@ /*=================================================================== 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 "mitkPAVolume.h" #include #include #include mitk::pa::Volume::Volume(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing) { + MITK_INFO << "Initialized by data*"; if (data == nullptr) mitkThrow() << "You may not initialize a mitk::Volume with a nullptr"; m_InternalMitkImage = mitk::Image::New(); unsigned int* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; dimensions[0] = yDim; dimensions[1] = xDim; dimensions[2] = zDim; m_XDim = xDim; m_YDim = yDim; m_ZDim = zDim; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); m_InternalMitkImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); m_InternalMitkImage->SetImportVolume(data, Image::ImportMemoryManagementType::CopyMemory); SetSpacing(spacing); m_FastAccessDataPointer = GetData(); delete data; } mitk::pa::Volume::Volume(mitk::Image::Pointer image) { + MITK_INFO << "Initialized by mitk::Image"; + if (image.IsNull()) mitkThrow() << "You may not initialize a mitk::Volume with a null reference to an mitk image"; unsigned int* dimensions = image->GetDimensions(); - m_XDim = dimensions[1]; m_YDim = dimensions[0]; + m_XDim = dimensions[1]; m_ZDim = dimensions[2]; m_InternalMitkImage = image; m_FastAccessDataPointer = GetData(); } double mitk::pa::Volume::GetSpacing() { return m_InternalMitkImage->GetGeometry()->GetSpacing()[0]; } void mitk::pa::Volume::SetSpacing(double spacing) { const mitk::ScalarType spacingArray[]{ spacing, spacing, spacing }; m_InternalMitkImage->SetSpacing(spacingArray); } mitk::pa::Volume::~Volume() { m_InternalMitkImage = nullptr; } mitk::pa::Volume::Pointer mitk::pa::Volume::New(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing) { mitk::pa::Volume::Pointer smartPtr = new mitk::pa::Volume(data, xDim, yDim, zDim, spacing); smartPtr->UnRegister(); return smartPtr; } mitk::pa::Volume::Pointer mitk::pa::Volume::New(mitk::Image::Pointer image) { mitk::pa::Volume::Pointer smartPtr = new mitk::pa::Volume(image); smartPtr->UnRegister(); return smartPtr; } mitk::Image::Pointer mitk::pa::Volume::AsMitkImage() { return m_InternalMitkImage; } mitk::pa::Volume::Pointer mitk::pa::Volume::DeepCopy() { long length = GetXDim()*GetYDim()*GetZDim(); double* data = new double[length]; memcpy(data, GetData(), length * sizeof(double)); return mitk::pa::Volume::New(data, GetXDim(), GetYDim(), GetZDim(), GetSpacing()); } double mitk::pa::Volume::GetData(unsigned int x, unsigned int y, unsigned int z) { return m_FastAccessDataPointer[GetIndex(x, y, z)]; } void mitk::pa::Volume::SetData(double data, unsigned int x, unsigned int y, unsigned int z) { m_FastAccessDataPointer[GetIndex(x, y, z)] = data; } unsigned int mitk::pa::Volume::GetXDim() { return m_XDim; } unsigned int mitk::pa::Volume::GetYDim() { return m_YDim; } unsigned int mitk::pa::Volume::GetZDim() { return m_ZDim; } double* mitk::pa::Volume::GetData() const { mitk::ImageWriteAccessor imgRead(m_InternalMitkImage, m_InternalMitkImage->GetVolumeData()); return (double*)imgRead.GetData(); } int mitk::pa::Volume::GetIndex(unsigned int x, unsigned int y, unsigned int z) { #ifdef _DEBUG if (x < 0 || x >(GetXDim() - 1) || y < 0 || y >(GetYDim() - 1) || z < 0 || z >(GetZDim() - 1)) { MITK_ERROR << "Index out of bounds at " << x << "|" << y << "|" << z; mitkThrow() << "Index out of bounds exception!"; } #endif return z * m_XDim * m_YDim + x * m_YDim + y; } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp index f4f49f6d4b..a0193346eb 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp @@ -1,98 +1,85 @@ /*=================================================================== 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 "mitkPAVesselDrawer.h" #include "mitkPAVesselProperties.h" mitk::pa::VesselDrawer::VesselDrawer() { } mitk::pa::VesselDrawer::~VesselDrawer() { } void mitk::pa::VesselDrawer::DrawVesselInVolume( VesselProperties::Pointer properties, InSilicoTissueVolume::Pointer volume) { Vector::Pointer stepDirection = properties->GetDirectionVector(); Vector::Pointer fromPosition = properties->GetPositionVector()->Clone(); Vector::Pointer toPosition = properties->GetPositionVector()->Clone(); Vector::Pointer totalWalkingDistance = stepDirection->Clone(); totalWalkingDistance->Scale(1.0 / volume->GetSpacing()); MITK_INFO << "STEP DIRECTION NORM: " << stepDirection->GetNorm(); MITK_INFO << "TOTAL WALKING DISTANCE NORM: " << totalWalkingDistance->GetNorm(); while (totalWalkingDistance->GetNorm() >= 1) { double xPos = fromPosition->GetElement(0); double yPos = fromPosition->GetElement(1); double zPos = fromPosition->GetElement(2); if (!volume->IsInsideVolume(xPos, yPos, zPos)) { properties->SetRadiusInVoxel(0); return; } double radius = properties->GetRadiusInVoxel(); double ceiledRadius = ceil(radius); for (int x = xPos - ceiledRadius; x <= xPos + ceiledRadius; x += 1) for (int y = yPos - ceiledRadius; y <= yPos + ceiledRadius; y += 1) for (int z = zPos - ceiledRadius; z <= zPos + ceiledRadius; z += 1) { if (!volume->IsInsideVolume(x, y, z)) { continue; } double xDiff = x - xPos; double yDiff = y - yPos; double zDiff = z - zPos; double vectorLengthDiff = radius - sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); if (vectorLengthDiff > 0) { volume->SetVolumeValues(x, y, z, properties->GetAbsorptionCoefficient(), properties->GetScatteringCoefficient(), properties->GetAnisotopyCoefficient(), mitk::pa::InSilicoTissueVolume::SegmentationType::VESSEL); } } totalWalkingDistance->Subtract(stepDirection); fromPosition->Add(stepDirection); } properties->SetPositionVector(fromPosition); } - -bool mitk::pa::Equal(const VesselDrawer::Pointer leftHandSide, const VesselDrawer::Pointer rightHandSide, double eps, bool verbose) -{ - MITK_INFO(verbose) << "=== mitk::pa::Vessel Equal ==="; - - if (rightHandSide.IsNull() || leftHandSide.IsNull()) - { - MITK_INFO(verbose) << "Cannot compare nullpointers"; - return false; - } - - return true; -} diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp index 153a15be90..9a2abb9380 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp @@ -1,220 +1,225 @@ /*=================================================================== 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 "mitkPAVolumeManipulator.h" #include "mitkPAExceptions.h" #include "mitkPAInSilicoTissueVolume.h" #include "itkResampleImageFilter.h" +#include "itkGaussianInterpolateImageFunction.h" // Includes for image casting between ITK and MITK #include #include +#include + mitk::pa::VolumeManipulator::VolumeManipulator() {} mitk::pa::VolumeManipulator::~VolumeManipulator() {} void mitk::pa::VolumeManipulator::ThresholdImage(Volume::Pointer image, double threshold) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) if (image->GetData(x, y, z) > threshold) image->SetData(1, x, y, z); else image->SetData(0, x, y, z); } void mitk::pa::VolumeManipulator::MultiplyImage(Volume::Pointer image, double factor) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) image->SetData(image->GetData(x, y, z)*factor, x, y, z); } void mitk::pa::VolumeManipulator::Log10Image(Volume::Pointer image) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) { if (image->GetData(x, y, z) < 0) { MITK_ERROR << "Signal was smaller than 0. There is an error in the input image file."; throw InvalidValueException("Signal was smaller than 0. There is an error in the input image file."); } image->SetData(log10(image->GetData(x, y, z)), x, y, z); } } void mitk::pa::VolumeManipulator::RescaleImage(InSilicoTissueVolume::Pointer image, double ratio) { - image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio)); image->SetScatteringVolume(RescaleImage(image->GetScatteringVolume(), ratio)); image->SetAnisotropyVolume(RescaleImage(image->GetAnisotropyVolume(), ratio)); image->SetSegmentationVolume(RescaleImage(image->GetSegmentationVolume(), ratio)); + image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio)); } mitk::pa::Volume::Pointer mitk::pa::VolumeManipulator::RescaleImage(Volume::Pointer image, double ratio) { typedef itk::Image ImageType; typedef itk::ResampleImageFilter FilterType; + typedef itk::GaussianInterpolateImageFunction InterpolatorType; auto input = image->AsMitkImage(); ImageType::Pointer itkInput = ImageType::New(); mitk::CastToItkImage(input, itkInput); ImageType::SizeType outputSize; outputSize[0] = input->GetDimensions()[0] * ratio; outputSize[1] = input->GetDimensions()[1] * ratio; outputSize[2] = input->GetDimensions()[2] * ratio; FilterType::Pointer resampleImageFilter = FilterType::New(); resampleImageFilter->SetInput(itkInput); resampleImageFilter->SetSize(outputSize); - resampleImageFilter->SetOutputSpacing(input->GetGeometry()->GetSpacing()[0] * ratio); + resampleImageFilter->SetInterpolator(InterpolatorType::New()); + resampleImageFilter->SetOutputSpacing(input->GetGeometry()->GetSpacing()[0] / ratio); resampleImageFilter->UpdateLargestPossibleRegion(); ImageType::Pointer output = resampleImageFilter->GetOutput(); mitk::Image::Pointer mitkOutput = mitk::Image::New(); GrabItkImageMemory(output, mitkOutput); return Volume::New(mitkOutput); } /** * @brief Fast 3D Gaussian convolution IIR approximation * @param paVolume * @param sigma * @author Pascal Getreuer * * Copyright (c) 2011, Pascal Getreuer * All rights reserved. * * This program is free software: you can redistribute it and/or modify it * under the terms of the simplified BSD license. * * You should have received a copy of these licenses along with this program. * If not, see . */ void mitk::pa::VolumeManipulator::GaussianBlur3D(mitk::pa::Volume::Pointer paVolume, double sigma) { double* volume = paVolume->GetData(); long width = paVolume->GetYDim(); long height = paVolume->GetXDim(); long depth = paVolume->GetZDim(); const long plane = width*height; const long numel = plane*depth; double lambda, dnu; double nu, boundaryscale, postscale; double *ptr; long i, x, y, z; int step; if (sigma <= 0) return; lambda = (sigma*sigma) / (8.0); dnu = (1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda)) / (2.0*lambda); nu = dnu; boundaryscale = 1.0 / (1.0 - dnu); postscale = pow(dnu / lambda, 12); /* Filter horizontally along each row */ for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { for (step = 0; step < 4; step++) { ptr = volume + width*(y + height*z); ptr[0] *= boundaryscale; /* Filter rightwards */ for (x = 1; x < width; x++) { ptr[x] += nu*ptr[x - 1]; } ptr[x = width - 1] *= boundaryscale; /* Filter leftwards */ for (; x > 0; x--) { ptr[x - 1] += nu*ptr[x]; } } } } /* Filter vertically along each column */ for (z = 0; z < depth; z++) { for (x = 0; x < width; x++) { for (step = 0; step < 4; step++) { ptr = volume + x + plane*z; ptr[0] *= boundaryscale; /* Filter downwards */ for (i = width; i < plane; i += width) { ptr[i] += nu*ptr[i - width]; } ptr[i = plane - width] *= boundaryscale; /* Filter upwards */ for (; i > 0; i -= width) { ptr[i - width] += nu*ptr[i]; } } } } /* Filter along z-dimension */ for (y = 0; y < height; y++) { for (x = 0; x < width; x++) { for (step = 0; step < 4; step++) { ptr = volume + x + width*y; ptr[0] *= boundaryscale; for (i = plane; i < numel; i += plane) { ptr[i] += nu*ptr[i - plane]; } ptr[i = numel - plane] *= boundaryscale; for (; i > 0; i -= plane) { ptr[i - plane] += nu*ptr[i]; } } } } for (i = 0; i < numel; i++) { volume[i] *= postscale; } }