diff --git a/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h b/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h index 3f1b242720..2b3575fbb1 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h +++ b/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h @@ -1,157 +1,160 @@ /*=================================================================== 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 MITKPHOTOACOUSTICVOLUME_H #define MITKPHOTOACOUSTICVOLUME_H #include #include #include #include #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT InSilicoTissueVolume : public itk::LightObject { public: mitkClassMacroItkParent(InSilicoTissueVolume, itk::LightObject) mitkNewMacro1Param(Self, TissueGeneratorParameters::Pointer) enum SegmentationType { AIR = -1, BACKGROUND = 0, VESSEL = 1, FAT = 2, SKIN = 3 }; /** * @brief ConvertToMitkImage * @return a pointer to an mitk image containing this volume. */ mitk::Image::Pointer ConvertToMitkImage(); /** * @brief SetVolumeValues sets the values for aborption, scattering and anisotropy at the specified voxel location. * * @param x * @param y * @param z * @param absorption * @param scattering * @param anisotropy * @param segmentType */ void SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType); /** * @brief SetVolumeValues sets the values for aborption, scattering and anisotropy at the specified voxel location. * * @param x * @param y * @param z * @param absorption * @param scattering * @param anisotropy */ void SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy); /** * @brief IsInsideVolume * * @param x * @param y * @param z * @return true if the voxel location is inside the volume */ bool IsInsideVolume(int x, int y, int z); /** * @brief AddDoubleProperty adds a persistent property to the volume, which will be exported to the mitk image. * * @param label * @param value */ void AddDoubleProperty(std::string label, double value); /** * @brief AddIntProperty adds a persistent property to the volume, which will be exported to the mitk image. * * @param label * @param value */ void AddIntProperty(std::string label, int value); Volume::Pointer GetAbsorptionVolume(); Volume::Pointer GetScatteringVolume(); Volume::Pointer GetAnisotropyVolume(); Volume::Pointer GetSegmentationVolume(); void SetAbsorptionVolume(Volume::Pointer volume); void SetScatteringVolume(Volume::Pointer volume); void SetAnisotropyVolume(Volume::Pointer volume); void SetSegmentationVolume(Volume::Pointer volume); + double GetSpacing(); + void SetSpacing(double spacing); + void FinalizeVolume(); itkGetMacro(TissueParameters, TissueGeneratorParameters::Pointer); itkGetMacro(TDim, unsigned int); static InSilicoTissueVolume::Pointer New(mitk::pa::Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList); protected: InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters); InSilicoTissueVolume(Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList); virtual ~InSilicoTissueVolume(); mitk::pa::Volume::Pointer m_AbsorptionVolume; mitk::pa::Volume::Pointer m_ScatteringVolume; mitk::pa::Volume::Pointer m_AnisotropyVolume; mitk::pa::Volume::Pointer m_SegmentationVolume; TissueGeneratorParameters::Pointer m_TissueParameters; unsigned int m_TDim; void RandomizeTissueCoefficients(long rngSeed, bool useRngSeed, double percentage); mitk::PropertyList::Pointer m_PropertyList; private: void FillZLayer(int x, int y, double startIdx, double endIdx, double absorption, double scattering, double anisotropy, SegmentationType segmentationType); void AddSkinAndAirLayers(); void UpdatePropertyList(); }; } } #endif // MITKPHOTOACOUSTICVOLUME_H diff --git a/Modules/PhotoacousticsLib/include/mitkPAVector.h b/Modules/PhotoacousticsLib/include/mitkPAVector.h index 4ef09eadef..0fb74e31a3 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVector.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVector.h @@ -1,132 +1,136 @@ /*=================================================================== 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 MITKSMARTVECTOR_H #define MITKSMARTVECTOR_H #include #include #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT Vector : public itk::LightObject { public: mitkClassMacroItkParent(Vector, itk::LightObject) itkFactorylessNewMacro(Self) /** * @brief GetNorm calculates the length of this vector. * @return the euclidean norm */ double GetNorm(); double GetElement(unsigned short index); void SetElement(unsigned short index, double value); /** * @brief Normalize normalizes this vector. After calling this GetNorm() will return 1. */ void Normalize(); void SetValue(Vector::Pointer value); /** * @brief RandomizeByPercentage alters this vector randomly by [-percentage, percentage] of the bendingFactor. * * @param percentage * @param bendingFactor */ void RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937* rng); /** * @brief Randomize randomizes this vector to be [lowerLimit, upperLimit] in each element * * @param xLowerLimit * @param xUpperLimit * @param yLowerLimit * @param yUpperLimit * @param zLowerLimit * @param zUpperLimit */ void Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937* rng); /** * @brief Randomize randomizes this vector to be [0, limit] in each element * * @param xLimit * @param yLimit * @param zLimit */ void Randomize(double xLimit, double yLimit, double zLimit, std::mt19937* rng); /** * @brief Randomize randomizes this vector to be [-1, 1] in each element */ void Randomize(std::mt19937* rng); /** * @brief Rotate rotates this Vector around the x, y and z axis with the given angles in radians * * @param thetaChange rotation of the inclination angle in radians * @param phiChange rotation of the azimuthal angle in radians */ void Rotate(double xAngle, double yAngle); /** * @brief Scale scales this Vector with the given factor * * @param factor the scaling factor * * If a negative number is provided, the direction of the vector will be inverted. */ void Scale(double factor); /** * @brief Clone create a deep copy of this vector * * @return a new vector with the same values. */ Vector::Pointer Clone(); + void Subtract(Vector::Pointer other); + + void Add(Vector::Pointer other); + protected: Vector(); virtual ~Vector(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const override; private: mitk::Vector3D m_Vector; }; /** * @brief Equal A function comparing two vectors for beeing equal * * @param rightHandSide A Vector to be compared * @param leftHandSide A Vector 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 Vector::Pointer leftHandSide, const Vector::Pointer rightHandSide, double eps, bool verbose); } } #endif // MITKSMARTVECTOR_H diff --git a/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h index f82ccf104b..7b1332a8c7 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h @@ -1,63 +1,63 @@ /*=================================================================== 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(Vector::Pointer toPosition, + 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/include/mitkPAVolume.h b/Modules/PhotoacousticsLib/include/mitkPAVolume.h index bfbad861c3..9a6f9c9b59 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVolume.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVolume.h @@ -1,145 +1,149 @@ /*=================================================================== 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 MITKPHOTOACOUSTIC3dVOLUME_H #define MITKPHOTOACOUSTIC3dVOLUME_H #include "MitkPhotoacousticsLibExports.h" //Includes for smart pointer usage #include #include namespace mitk { namespace pa { /** * @brief The Volume class is designed to encapsulate volumetric information and to provide convenience methods * for data access and image conversions. */ class MITKPHOTOACOUSTICSLIB_EXPORT Volume : public itk::LightObject { public: mitkClassMacroItkParent(Volume, itk::LightObject); /** *@brief returns smartpointer reference to a new instance of this objects. * The given data array will be freed upon calling this constructor. *@param data *@param xDim *@param yDim *@param zDim *@return smartpointer reference to a new instance of this object */ - static Volume::Pointer New(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim); + static Volume::Pointer New(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing); static Volume::Pointer New(mitk::Image::Pointer image); /** * @brief GetData. Returns data at wanted position. For performance reasons, this method will not check, * if the specified position it within the array. Please use the GetXDim(), GetYDim() and GetZDim() methods * to check for yourself if necessary. * * @param x * @param y * @param z * @return the data contained within the data array held by this Volume at * positon x|y|z. */ double GetData(unsigned int x, unsigned int y, unsigned int z); /** * Returns a const reference to the data encapsuled by this class. */ double* GetData() const; /** * @brief SetData * @param data * @param x * @param y * @param z */ void SetData(double data, unsigned int x, unsigned int y, unsigned int z); /** * @brief GetXDim * @return size of x dimension of this Volume */ unsigned int GetXDim(); /** * @brief GetYDim * @return size of y dimension of this Volume */ unsigned int GetYDim(); /** * @brief GetZDim * @return size of z dimension of this Volume */ unsigned int GetZDim(); /** *@brief returns the Volume instance as an mitk image */ Image::Pointer AsMitkImage(); /** * @brief DeepCopy * @return a deep copy of this Volume. the old volume remains intact and memory is NOT shared * between the objects. */ Volume::Pointer DeepCopy(); /** *@brief convenience method to enable consistent access to the dat array *@return a 1d index from 3d pixel coordinates */ int GetIndex(unsigned int x, unsigned int y, unsigned int z); + double GetSpacing(); + + void SetSpacing(double spacing); + protected: /** * @brief Initialize initializes this volume with the given pointer to the data array. * It is assumed, that the array is of dimension xDim|yDim|zDim. * The Photoacoustic3DVolume will handle memory management of the array and delete it on * constructor call. * * @param data a pointer to the data array * @param xDim x dimension of the data * @param yDim y dimension of the data * @param zDim z dimension of the data */ - Volume(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim); + Volume(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing); Volume(mitk::Image::Pointer image); virtual ~Volume(); const int NUMBER_OF_SPATIAL_DIMENSIONS = 3; Image::Pointer m_InternalMitkImage; // this data is kept to enable fast access unsigned int m_XDim; unsigned int m_YDim; unsigned int m_ZDim; double* m_FastAccessDataPointer; }; } } #endif // MITKPHOTOACOUSTIC3dVOLUME_H diff --git a/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h b/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h index dd42d19565..4b1eada957 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h @@ -1,58 +1,58 @@ /*=================================================================== 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 MITKPHOTOACOUSTIC3DVOLUMEMANIPULATOR_H #define MITKPHOTOACOUSTIC3DVOLUMEMANIPULATOR_H #include #include #include "mitkPAVolume.h" #include "mitkPAInSilicoTissueVolume.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT VolumeManipulator final { public: /** * @brief ThresholdImage applies a binary threshold filter to this image. * @param threshold */ - static void ThresholdImage(mitk::pa::Volume::Pointer image, double threshold); + static void ThresholdImage(Volume::Pointer image, double threshold); /** * @brief Multiplies the image with a given factor * @param factor */ - static void MultiplyImage(mitk::pa::Volume::Pointer image, double factor); + static void MultiplyImage(Volume::Pointer image, double factor); - static void GaussianBlur3D(mitk::pa::Volume::Pointer paVolume, double sigma); + static void GaussianBlur3D(Volume::Pointer paVolume, double sigma); - static void Log10Image(mitk::pa::Volume::Pointer image); + static void Log10Image(Volume::Pointer image); static void RescaleImage(InSilicoTissueVolume::Pointer image, double ratio); - static mitk::pa::Volume::Pointer RescaleImage(Volume::Pointer image, double ratio); + static Volume::Pointer RescaleImage(Volume::Pointer image, double ratio); private: VolumeManipulator(); virtual ~VolumeManipulator(); }; } } #endif // MITKPHOTOACOUSTIC3DVOLUMEMANIPULATOR_H diff --git a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp index 09cc8dc1f3..7e9090acb7 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp @@ -1,118 +1,117 @@ /*=================================================================== 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 "mitkPAVessel.h" #include #define _USE_MATH_DEFINES #include #include mitk::pa::Vessel::Vessel(VesselProperties::Pointer initialProperties) : m_RangeDistribution(M_PI / 16, M_PI / 8), m_SignDistribution(-1, 1) { m_Finished = false; //Copy this so it may be reused for other vessels. m_VesselProperties = VesselProperties::New(initialProperties); m_RadiusRangeDistribution = std::uniform_real_distribution<>(NEW_RADIUS_MINIMUM_RELATIVE_SIZE, NEW_RADIUS_MAXIMUM_RELATIVE_SIZE); m_VesselMeanderStrategy = VesselMeanderStrategy::New(); m_WalkedDistance = 0; m_VesselDrawer = VesselDrawer::New(); } mitk::pa::Vessel::~Vessel() { m_VesselProperties = nullptr; m_VesselMeanderStrategy = nullptr; } void mitk::pa::Vessel::ExpandVessel(InSilicoTissueVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng) { - Vector::Pointer oldPosition = m_VesselProperties->GetPositionVector()->Clone(); + m_VesselDrawer->DrawVesselInVolume(m_VesselProperties, volume); (m_VesselMeanderStrategy->*calculateNewPosition)(m_VesselProperties->GetPositionVector(), m_VesselProperties->GetDirectionVector(), bendingFactor, rng); - m_VesselDrawer->DrawVesselInVolume(oldPosition, m_VesselProperties, volume); m_WalkedDistance += m_VesselProperties->GetDirectionVector()->GetNorm(); } bool mitk::pa::Vessel::CanBifurcate() { return m_VesselProperties->GetBifurcationFrequency() < m_WalkedDistance; } int mitk::pa::Vessel::GetSign(std::mt19937 *rng) { if (m_SignDistribution(*rng) < 0) return -1; return 1; } mitk::pa::Vessel::Pointer mitk::pa::Vessel::Bifurcate(std::mt19937* rng) { VesselProperties::Pointer vesselParams = VesselProperties::New(m_VesselProperties); double thetaChange = m_RangeDistribution(*rng) * GetSign(rng); double phiChange = m_RangeDistribution(*rng) * GetSign(rng); vesselParams->GetDirectionVector()->Rotate(thetaChange, phiChange); m_VesselProperties->GetDirectionVector()->Rotate(-thetaChange, -phiChange); double newRadius = m_RadiusRangeDistribution(*rng)*m_VesselProperties->GetRadiusInVoxel(); vesselParams->SetRadiusInVoxel(newRadius); m_VesselProperties->SetRadiusInVoxel( sqrt(m_VesselProperties->GetRadiusInVoxel()*m_VesselProperties->GetRadiusInVoxel() - newRadius*newRadius)); m_WalkedDistance = 0; return Vessel::New(vesselParams); } bool mitk::pa::Vessel::IsFinished() { return m_VesselProperties->GetRadiusInVoxel() < MINIMUM_VESSEL_RADIUS; } bool mitk::pa::Equal(const Vessel::Pointer leftHandSide, const Vessel::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; } if (leftHandSide->IsFinished() != rightHandSide->IsFinished()) { MITK_INFO(verbose) << "Not same finished state."; return false; } if (leftHandSide->CanBifurcate() != rightHandSide->CanBifurcate()) { MITK_INFO(verbose) << "Not same bifurcation state."; return false; } if (!Equal(leftHandSide->GetVesselProperties(), rightHandSide->GetVesselProperties(), eps, verbose)) { MITK_INFO(verbose) << "Vesselproperties not equal"; return false; } return true; } diff --git a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp index 679f0bda4a..6ddf981e18 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp @@ -1,59 +1,47 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPAVesselMeanderStrategy.h" mitk::pa::VesselMeanderStrategy::VesselMeanderStrategy() { } mitk::pa::VesselMeanderStrategy::~VesselMeanderStrategy() { } void mitk::pa::VesselMeanderStrategy::CalculateNewPositionInStraightLine( Vector::Pointer position, Vector::Pointer direction, double /*bendingFactor*/, std::mt19937* rng) { if (direction->GetNorm() <= mitk::eps) { direction->Randomize(rng); } - - direction->Normalize(); - - position->SetElement(0, position->GetElement(0) + direction->GetElement(0)); - position->SetElement(1, position->GetElement(1) + direction->GetElement(1)); - position->SetElement(2, position->GetElement(2) + direction->GetElement(2)); - - position->Normalize(); } void mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition( Vector::Pointer position, Vector::Pointer direction, double bendingFactor, std::mt19937* rng) { if (direction->GetNorm() <= mitk::eps) { direction->Randomize(rng); } direction->RandomizeByPercentage(RANDOMIZATION_PERCENTAGE, bendingFactor, rng); direction->Normalize(); - - position->SetElement(0, position->GetElement(0) + direction->GetElement(0)); - position->SetElement(1, position->GetElement(1) + direction->GetElement(1)); - position->SetElement(2, position->GetElement(2) + direction->GetElement(2)); } diff --git a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp index 2666b3c3c1..9c9e38f29b 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp @@ -1,352 +1,365 @@ /*=================================================================== 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); - m_ScatteringVolume = Volume::New(scatteringArray, xDim, yDim, zDim); - m_AnisotropyVolume = Volume::New(anisotropyArray, xDim, yDim, zDim); - m_SegmentationVolume = Volume::New(segmentationArray, xDim, yDim, zDim); + 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 70fb4f4944..fd65f072c4 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp @@ -1,137 +1,150 @@ /*=================================================================== 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) + unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing) { 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) { 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_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) +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); + 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()); + 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/Generator/mitkPASlicedVolumeGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp index ce6ea25487..b514c7f889 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp @@ -1,118 +1,118 @@ /*=================================================================== 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 "mitkPASlicedVolumeGenerator.h" #include mitk::pa::SlicedVolumeGenerator::SlicedVolumeGenerator(int centralYSlice, bool precorrect, Volume::Pointer precorrectionVolume, bool inverse) { m_CentralYSlice = centralYSlice; m_Precorrect = precorrect; m_Inverse = inverse; m_PrecorrectionVolume = nullptr; if (m_Precorrect) { m_PrecorrectionVolume = precorrectionVolume; } } mitk::pa::SlicedVolumeGenerator::~SlicedVolumeGenerator() { m_CentralYSlice = -1; m_Precorrect = false; m_PrecorrectionVolume = nullptr; } mitk::pa::Volume::Pointer mitk::pa::SlicedVolumeGenerator::GetSlicedFluenceImageFromComposedVolume( ComposedVolume::Pointer composedVolume) { int fluenceComponents = composedVolume->GetNumberOfFluenceComponents(); int xDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetXDim(); int zDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetZDim(); double* imageArray = new double[xDim*zDim*fluenceComponents]; for (int fluenceComponentIdx = 0; fluenceComponentIdx < fluenceComponents; fluenceComponentIdx++) for (int z = 0; z < zDim; z++) for (int x = 0; x < xDim; x++) { int index = z * xDim * fluenceComponents + x * fluenceComponents + fluenceComponentIdx; imageArray[index] = composedVolume->GetFluenceValue(fluenceComponentIdx, x, m_CentralYSlice, z); if (m_Precorrect) { imageArray[index] = imageArray[index] / m_PrecorrectionVolume->GetData(x, m_CentralYSlice, z); } if (m_Inverse) { if (abs(imageArray[index] - 0) >= mitk::eps) imageArray[index] = 1 / imageArray[index]; else imageArray[index] = INFINITY; } } - return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim); + return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim, composedVolume->GetGroundTruthVolume()->GetSpacing()); } mitk::pa::Volume::Pointer mitk::pa::SlicedVolumeGenerator::GetSlicedSignalImageFromComposedVolume( ComposedVolume::Pointer composedVolume) { int fluenceComponents = composedVolume->GetNumberOfFluenceComponents(); int xDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetXDim(); int zDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetZDim(); double* imageArray = new double[xDim*zDim*fluenceComponents]; for (int fluenceComponentIdx = 0; fluenceComponentIdx < fluenceComponents; fluenceComponentIdx++) for (int z = 0; z < zDim; z++) for (int x = 0; x < xDim; x++) { int y = m_CentralYSlice + composedVolume->GetYOffsetForFluenceComponentInPixels(fluenceComponentIdx); imageArray[z * xDim * fluenceComponents + x * fluenceComponents + fluenceComponentIdx] = composedVolume->GetFluenceValue(fluenceComponentIdx, x, m_CentralYSlice, z) * composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetData(x, y, z); } - return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim); + return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim, composedVolume->GetGroundTruthVolume()->GetSpacing()); } mitk::pa::Volume::Pointer mitk::pa::SlicedVolumeGenerator::GetSlicedGroundTruthImageFromComposedVolume( ComposedVolume::Pointer composedVolume) { int fluenceComponents = composedVolume->GetNumberOfFluenceComponents(); int xDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetXDim(); int zDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetZDim(); double* imageArray = new double[xDim*zDim*fluenceComponents]; for (int fluenceComponentIdx = 0; fluenceComponentIdx < fluenceComponents; fluenceComponentIdx++) for (int z = 0; z < zDim; z++) for (int x = 0; x < xDim; x++) { int y = m_CentralYSlice + composedVolume->GetYOffsetForFluenceComponentInPixels(fluenceComponentIdx); imageArray[z * xDim * fluenceComponents + x * fluenceComponents + fluenceComponentIdx] = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetData(x, y, z); } - return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim); + return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim, composedVolume->GetGroundTruthVolume()->GetSpacing()); } diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp index a3ecd9c092..41ce37435e 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp @@ -1,176 +1,176 @@ /*=================================================================== 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 "mitkPATissueGenerator.h" #include "mitkPAVector.h" #include "mitkPAVolumeManipulator.h" mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData( TissueGeneratorParameters::Pointer parameters) { MITK_DEBUG << "Initializing Empty Volume"; + if (parameters->GetDoPartialVolume()) + { + parameters->SetXDim(parameters->GetXDim() * 10); + parameters->SetYDim(parameters->GetYDim() * 10); + parameters->SetZDim(parameters->GetZDim() * 10); + parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() * 10); + parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() / 10); + } + auto generatedVolume = mitk::pa::InSilicoTissueVolume::New(parameters); const double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; std::mt19937 randomNumberGenerator; std::random_device randomDevice; if (parameters->GetUseRngSeed()) { randomNumberGenerator.seed(parameters->GetRngSeed()); } else { if (randomDevice.entropy() > 0.1) { randomNumberGenerator.seed(randomDevice()); } else { randomNumberGenerator.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); } } - if (parameters->GetDoPartialVolume()) - { - parameters->SetXDim(parameters->GetXDim() * 10); - parameters->SetYDim(parameters->GetYDim() * 10); - parameters->SetZDim(parameters->GetZDim() * 10); - parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() * 10); - parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() / 10); - } - std::uniform_int_distribution randomNumVesselDistribution(parameters->GetMinNumberOfVessels(), parameters->GetMaxNumberOfVessels()); std::uniform_real_distribution randomBendingDistribution(parameters->GetMinVesselBending(), parameters->GetMaxVesselBending()); std::uniform_real_distribution randomAbsorptionDistribution(parameters->GetMinVesselAbsorption(), parameters->GetMaxVesselAbsorption()); std::uniform_real_distribution randomRadiusDistribution(parameters->GetMinVesselRadiusInMillimeters(), parameters->GetMaxVesselRadiusInMillimeters()); std::uniform_real_distribution randomScatteringDistribution(parameters->GetMinVesselScattering(), parameters->GetMaxVesselScattering()); std::uniform_real_distribution randomAnisotropyDistribution(parameters->GetMinVesselAnisotropy(), parameters->GetMaxVesselAnisotropy()); std::uniform_int_distribution borderTypeDistribution(0, 3); int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator); generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels); generatedVolume->AddIntProperty("bifurcationFrequency", parameters->GetVesselBifurcationFrequency()); MITK_INFO << "Simulating " << numberOfBloodVessels << " vessel structures"; for (int vesselNumber = 0; vesselNumber < numberOfBloodVessels; vesselNumber++) { Vector::Pointer initialPosition = Vector::New(); Vector::Pointer initialDirection = Vector::New(); double initialRadius = randomRadiusDistribution(randomNumberGenerator) / parameters->GetVoxelSpacingInCentimeters() / 10; std::stringstream radiusString; radiusString << "vessel_" << vesselNumber + 1 << "_radius"; generatedVolume->AddDoubleProperty(radiusString.str(), initialRadius); double absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator); std::stringstream absorptionString; absorptionString << "vessel_" << vesselNumber + 1 << "_absorption"; generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient); double bendingFactor = randomBendingDistribution(randomNumberGenerator); std::stringstream bendingString; bendingString << "vessel_" << vesselNumber + 1 << "_bendingFactor"; generatedVolume->AddDoubleProperty(bendingString.str(), bendingFactor); double vesselScattering = randomScatteringDistribution(randomNumberGenerator); std::stringstream scatteringString; scatteringString << "vessel_" << vesselNumber + 1 << "_scattering"; generatedVolume->AddDoubleProperty(scatteringString.str(), vesselScattering); double vesselAnisotropy = randomAnisotropyDistribution(randomNumberGenerator); std::stringstream anisotropyString; anisotropyString << "vessel_" << vesselNumber + 1 << "_anisotropy"; generatedVolume->AddDoubleProperty(anisotropyString.str(), vesselAnisotropy); /*The vessel tree shall start at one of the 4 sides of the volume. * The vessels will always be completely contained in the volume * when starting to meander. * They will meander in a direction perpendicular to the * plane they started from (within the limits of the * DIRECTION_VECTOR_INITIAL_VARIANCE) */ int borderType = borderTypeDistribution(randomNumberGenerator); switch (borderType) { case 0: initialPosition->Randomize(0, 0, initialRadius, parameters->GetYDim() - initialRadius, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(1, 2, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; case 1: initialPosition->Randomize(parameters->GetXDim(), parameters->GetXDim(), initialRadius, parameters->GetYDim() - initialRadius, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; case 2: initialPosition->Randomize(initialRadius, parameters->GetXDim() - initialRadius, 0, 0, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, 1, 2, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; case 3: initialPosition->Randomize(initialRadius, parameters->GetXDim() - initialRadius, parameters->GetYDim(), parameters->GetYDim(), parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; } VesselProperties::Pointer vesselParams = VesselProperties::New(); vesselParams->SetDirectionVector(initialDirection); vesselParams->SetPositionVector(initialPosition); vesselParams->SetRadiusInVoxel(initialRadius); vesselParams->SetAbsorptionCoefficient(absorptionCoefficient); vesselParams->SetScatteringCoefficient(vesselScattering); vesselParams->SetAnisotopyCoefficient(vesselAnisotropy); vesselParams->SetBifurcationFrequency(parameters->GetVesselBifurcationFrequency()); vesselParams->SetDoPartialVolume(parameters->GetDoPartialVolume()); VesselTree::Pointer vesselTree = VesselTree::New(vesselParams); while (!vesselTree->IsFinished()) { vesselTree->Step(generatedVolume, parameters->GetCalculateNewVesselPositionCallback(), bendingFactor, &randomNumberGenerator); } } if (parameters->GetDoPartialVolume()) { VolumeManipulator::RescaleImage(generatedVolume, 0.1); parameters->SetXDim(parameters->GetXDim() / 10); parameters->SetYDim(parameters->GetYDim() / 10); parameters->SetZDim(parameters->GetZDim() / 10); parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() / 10); parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() * 10); } generatedVolume->FinalizeVolume(); return generatedVolume; } mitk::pa::InSilicoTissueGenerator::InSilicoTissueGenerator() { } mitk::pa::InSilicoTissueGenerator::~InSilicoTissueGenerator() { } diff --git a/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp b/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp index 3d63e78a68..b4d8709e64 100644 --- a/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp +++ b/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp @@ -1,258 +1,240 @@ #include "mitkPAIOUtil.h" #include "mitkIOUtil.h" #include "mitkImageReadAccessor.h" #include #include #include #include "mitkPAComposedVolume.h" #include "mitkPASlicedVolumeGenerator.h" #include "mitkPANoiseGenerator.h" #include "mitkPAVolumeManipulator.h" #include #include #include static std::vector splitString(const std::string &s, const char* delim) { std::vector elems; std::stringstream ss(s); std::string item; while (std::getline(ss, item, *delim)) { int numb; std::stringstream(item) >> numb; elems.push_back(numb); } return elems; } bool mitk::pa::IOUtil::DoesFileHaveEnding(std::string const &fullString, std::string const &ending) { if (fullString.length() == 0 || ending.length() == 0 || fullString.length() < ending.length()) return false; return (0 == fullString.compare(fullString.length() - ending.length(), ending.length(), ending)); } mitk::pa::IOUtil::IOUtil() {} mitk::pa::IOUtil::~IOUtil() {} mitk::pa::Volume::Pointer mitk::pa::IOUtil::LoadNrrd(std::string filename, double blur) { if (filename.empty() || filename == "") return nullptr; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(filename); if (inputImage.IsNull()) return nullptr; - int xDim = inputImage->GetDimensions()[1]; - int yDim = inputImage->GetDimensions()[0]; - int zDim = inputImage->GetDimensions()[2]; + auto returnImage = Volume::New(inputImage); - mitk::ImageReadAccessor readAccess(inputImage, inputImage->GetVolumeData(0)); - - double* dataArray = new double[xDim*yDim*zDim]; - const double* srcData = (const double*)readAccess.GetData(); - memcpy(dataArray, srcData, xDim*yDim*zDim * sizeof(double)); - - auto returnImage = mitk::pa::Volume::New(dataArray, xDim, yDim, zDim); - - mitk::pa::VolumeManipulator::GaussianBlur3D(returnImage, blur); + VolumeManipulator::GaussianBlur3D(returnImage, blur); return returnImage; } std::map mitk::pa::IOUtil::LoadFluenceContributionMaps(std::string foldername, double blur, int* progress, bool doLog10) { std::map resultMap; itk::Directory::Pointer directoryHandler = itk::Directory::New(); directoryHandler->Load(foldername.c_str()); for (unsigned int fileIndex = 0, numFiles = directoryHandler->GetNumberOfFiles(); fileIndex < numFiles; ++fileIndex) { std::string filename = std::string(directoryHandler->GetFile(fileIndex)); if (itksys::SystemTools::FileIsDirectory(filename)) continue; if (!DoesFileHaveEnding(filename, ".nrrd")) continue; size_t s = filename.find("_p"); size_t e = filename.find("Fluence", s); std::string sub = filename.substr(s + 2, e - s - 2); std::vector coords = splitString(sub, ","); if (coords.size() != 3) { MITK_ERROR << "Some of the data to read was corrupted or did not match the " << "naming pattern *_pN,N,NFluence*.nrrd"; mitkThrow() << "Some of the data to read was corrupted or did not match the" << " naming pattern *_pN,N,NFluence*.nrrd"; } else { MITK_DEBUG << "Extracted coords: " << coords[0] << "|" << coords[1] << "|" << coords[2] << " from string " << sub; Volume::Pointer nrrdFile = LoadNrrd(foldername + filename, blur); if (doLog10) VolumeManipulator::Log10Image(nrrdFile); resultMap[Position{ coords[0], coords[2] }] = nrrdFile; *progress = *progress + 1; } } return resultMap; } int mitk::pa::IOUtil::GetNumberOfNrrdFilesInDirectory(std::string directory) { return GetListOfAllNrrdFilesInDirectory(directory).size(); } std::vector mitk::pa::IOUtil::GetListOfAllNrrdFilesInDirectory(std::string directory, bool keepFileFormat) { std::vector filenames; itk::Directory::Pointer directoryHandler = itk::Directory::New(); directoryHandler->Load(directory.c_str()); for (unsigned int fileIndex = 0, numFiles = directoryHandler->GetNumberOfFiles(); fileIndex < numFiles; ++fileIndex) { std::string filename = std::string(directoryHandler->GetFile(fileIndex)); if (itksys::SystemTools::FileIsDirectory(filename)) continue; if (!DoesFileHaveEnding(filename, ".nrrd")) continue; if (keepFileFormat) { filenames.push_back(filename); } else { filenames.push_back(filename.substr(0, filename.size() - 5)); } } return filenames; } std::vector mitk::pa::IOUtil::GetAllChildfoldersFromFolder(std::string folderPath) { std::vector returnVector; itksys::Directory directoryHandler; directoryHandler.Load(folderPath.c_str()); for (unsigned int fileIndex = 0, numFiles = directoryHandler.GetNumberOfFiles(); fileIndex < numFiles; ++fileIndex) { std::string foldername = std::string(directoryHandler.GetFile(fileIndex)); std::string filename = folderPath + "/" + foldername; if (itksys::SystemTools::FileIsDirectory(filename)) { if (foldername != std::string(".") && foldername != std::string("..")) { MITK_INFO << filename; returnVector.push_back(filename); } continue; } //If there is a nrrd file in the directory we assume that a bottom level directory was chosen. if (DoesFileHaveEnding(filename, ".nrrd")) { returnVector.clear(); returnVector.push_back(folderPath); return returnVector; } } return returnVector; } mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::IOUtil::LoadInSilicoTissueVolumeFromNrrdFile(std::string nrrdFile) { MITK_INFO << "Initializing ComposedVolume by nrrd..."; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(nrrdFile); auto tissueParameters = TissueGeneratorParameters::New(); unsigned int xDim = inputImage->GetDimensions()[1]; unsigned int yDim = inputImage->GetDimensions()[0]; unsigned int zDim = inputImage->GetDimensions()[2]; tissueParameters->SetXDim(xDim); tissueParameters->SetYDim(yDim); tissueParameters->SetZDim(zDim); double xSpacing = inputImage->GetGeometry(0)->GetSpacing()[1]; double ySpacing = inputImage->GetGeometry(0)->GetSpacing()[0]; double zSpacing = inputImage->GetGeometry(0)->GetSpacing()[2]; if ((xSpacing - ySpacing) > mitk::eps || (xSpacing - zSpacing) > mitk::eps || (ySpacing - zSpacing) > mitk::eps) { throw mitk::Exception("Cannot handle unequal spacing."); } tissueParameters->SetVoxelSpacingInCentimeters(xSpacing); mitk::PropertyList::Pointer propertyList = inputImage->GetPropertyList(); mitk::ImageReadAccessor readAccess0(inputImage, inputImage->GetVolumeData(0)); double* m_AbsorptionArray = new double[xDim*yDim*zDim]; memcpy(m_AbsorptionArray, readAccess0.GetData(), xDim*yDim*zDim * sizeof(double)); - auto absorptionVolume = Volume::New(m_AbsorptionArray, xDim, yDim, zDim); + auto absorptionVolume = Volume::New(m_AbsorptionArray, xDim, yDim, zDim, xSpacing); mitk::ImageReadAccessor readAccess1(inputImage, inputImage->GetVolumeData(1)); double* m_ScatteringArray = new double[xDim*yDim*zDim]; memcpy(m_ScatteringArray, readAccess1.GetData(), xDim*yDim*zDim * sizeof(double)); - auto scatteringVolume = Volume::New(m_ScatteringArray, xDim, yDim, zDim); + auto scatteringVolume = Volume::New(m_ScatteringArray, xDim, yDim, zDim, xSpacing); mitk::ImageReadAccessor readAccess2(inputImage, inputImage->GetVolumeData(2)); double* m_AnisotropyArray = new double[xDim*yDim*zDim]; memcpy(m_AnisotropyArray, readAccess2.GetData(), xDim*yDim*zDim * sizeof(double)); - auto anisotropyVolume = Volume::New(m_AnisotropyArray, xDim, yDim, zDim); + auto anisotropyVolume = Volume::New(m_AnisotropyArray, xDim, yDim, zDim, xSpacing); Volume::Pointer segmentationVolume; if (inputImage->GetDimension() == 4) { mitk::ImageReadAccessor readAccess3(inputImage, inputImage->GetVolumeData(3)); double* m_SegmentationArray = new double[xDim*yDim*zDim]; memcpy(m_SegmentationArray, readAccess3.GetData(), xDim*yDim*zDim * sizeof(double)); - segmentationVolume = Volume::New(m_SegmentationArray, xDim, yDim, zDim); + segmentationVolume = Volume::New(m_SegmentationArray, xDim, yDim, zDim, xSpacing); } return mitk::pa::InSilicoTissueVolume::New(absorptionVolume, scatteringVolume, anisotropyVolume, segmentationVolume, tissueParameters, propertyList); } mitk::pa::FluenceYOffsetPair::Pointer mitk::pa::IOUtil::LoadFluenceSimulation(std::string fluenceSimulation) { MITK_INFO << "Adding slice..."; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(fluenceSimulation); - mitk::ImageReadAccessor readAccess0(inputImage, inputImage->GetVolumeData(0)); - - unsigned int xDim = inputImage->GetDimensions()[1]; - unsigned int yDim = inputImage->GetDimensions()[0]; - unsigned int zDim = inputImage->GetDimensions()[2]; - int size = xDim*yDim*zDim; - double* fluenceArray = new double[size]; - memcpy(fluenceArray, readAccess0.GetData(), size * sizeof(double)); auto yOffsetProperty = inputImage->GetProperty("y-offset"); if (yOffsetProperty.IsNull()) mitkThrow() << "No \"y-offset\" property found in fluence file!"; std::string yOff = yOffsetProperty->GetValueAsString(); MITK_INFO << "Reading y Offset: " << yOff; #ifdef __linux__ std::replace(yOff.begin(), yOff.end(), '.', ','); #endif // __linux__ double yOffset = std::stod(yOff); MITK_INFO << "Converted offset " << yOffset; - return mitk::pa::FluenceYOffsetPair::New(mitk::pa::Volume::New(fluenceArray, xDim, yDim, zDim), yOffset); + return FluenceYOffsetPair::New(Volume::New(inputImage), yOffset); } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp index 2d468a8062..449795b1b9 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp @@ -1,168 +1,182 @@ /*=================================================================== 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 "mitkPAVector.h" #include "chrono" #include "math.h" mitk::pa::Vector::Vector() { m_Vector.Fill(0); } mitk::pa::Vector::~Vector() { m_Vector.Fill(0); } double mitk::pa::Vector::GetNorm() { return m_Vector.GetNorm(); } double mitk::pa::Vector::GetElement(unsigned short index) { return m_Vector.GetElement(index); } void mitk::pa::Vector::SetElement(unsigned short index, double value) { m_Vector.SetElement(index, value); } void mitk::pa::Vector::Normalize() { double norm = m_Vector.GetNorm(); m_Vector.SetElement(0, m_Vector.GetElement(0) / norm); m_Vector.SetElement(1, m_Vector.GetElement(1) / norm); m_Vector.SetElement(2, m_Vector.GetElement(2) / norm); } void mitk::pa::Vector::SetValue(mitk::pa::Vector::Pointer value) { m_Vector.SetElement(0, value->GetElement(0)); m_Vector.SetElement(1, value->GetElement(1)); m_Vector.SetElement(2, value->GetElement(2)); } void mitk::pa::Vector::RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937* rng) { std::uniform_real_distribution<> range(-percentage, percentage); m_Vector.SetElement(0, m_Vector.GetElement(0) + (bendingFactor * range(*rng))); m_Vector.SetElement(1, m_Vector.GetElement(1) + (bendingFactor * range(*rng))); m_Vector.SetElement(2, m_Vector.GetElement(2) + (bendingFactor * range(*rng))); } void mitk::pa::Vector::Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937* rng) { std::uniform_real_distribution<> rangeX(xLowerLimit, xUpperLimit); std::uniform_real_distribution<> rangeY(yLowerLimit, yUpperLimit); std::uniform_real_distribution<> rangeZ(zLowerLimit, zUpperLimit); m_Vector.SetElement(0, rangeX(*rng)); m_Vector.SetElement(1, rangeY(*rng)); m_Vector.SetElement(2, rangeZ(*rng)); } void mitk::pa::Vector::Randomize(double xLimit, double yLimit, double zLimit, std::mt19937* rng) { Randomize(0, xLimit, 0, yLimit, 0, zLimit, rng); } void mitk::pa::Vector::Randomize(std::mt19937* rng) { Randomize(-1, 1, -1, 1, -1, 1, rng); } void mitk::pa::Vector::PrintSelf(std::ostream& os, itk::Indent /*indent*/) const { os << "X: " << m_Vector.GetElement(0) << std::endl; os << "Y: " << m_Vector.GetElement(1) << std::endl; os << "Z: " << m_Vector.GetElement(2) << std::endl; os << "Length: " << m_Vector.GetNorm() << std::endl; } void mitk::pa::Vector::Rotate(double thetaChange, double phiChange) { MITK_DEBUG << "Vector before rotation: (" << GetElement(0) << "|" << GetElement(1) << "|" << GetElement(2) << ")"; if (thetaChange == 0 && phiChange == 0) return; double x = GetElement(0); double y = GetElement(1); double z = GetElement(2); double r = sqrt(x*x + y*y + z*z); if (r == 0) return; double theta = acos(z / r); double phi = atan2(y, x); theta += thetaChange; phi += phiChange; SetElement(0, r * sin(theta) * cos(phi)); SetElement(1, r * sin(theta) * sin(phi)); SetElement(2, r * cos(theta)); MITK_DEBUG << "Vector after rotation: (" << GetElement(0) << "|" << GetElement(1) << "|" << GetElement(2) << ")"; } void mitk::pa::Vector::Scale(double factor) { m_Vector.SetElement(0, m_Vector.GetElement(0)*factor); m_Vector.SetElement(1, m_Vector.GetElement(1)*factor); m_Vector.SetElement(2, m_Vector.GetElement(2)*factor); } mitk::pa::Vector::Pointer mitk::pa::Vector::Clone() { auto returnVector = Vector::New(); returnVector->SetElement(0, this->GetElement(0)); returnVector->SetElement(1, this->GetElement(1)); returnVector->SetElement(2, this->GetElement(2)); return returnVector; } +void mitk::pa::Vector::Subtract(Vector::Pointer other) +{ + m_Vector.SetElement(0, m_Vector.GetElement(0) - other->GetElement(0)); + m_Vector.SetElement(1, m_Vector.GetElement(1) - other->GetElement(1)); + m_Vector.SetElement(2, m_Vector.GetElement(2) - other->GetElement(2)); +} + +void mitk::pa::Vector::Add(Vector::Pointer other) +{ + m_Vector.SetElement(0, m_Vector.GetElement(0) + other->GetElement(0)); + m_Vector.SetElement(1, m_Vector.GetElement(1) + other->GetElement(1)); + m_Vector.SetElement(2, m_Vector.GetElement(2) + other->GetElement(2)); +} + bool mitk::pa::Equal(const Vector::Pointer leftHandSide, const Vector::Pointer rightHandSide, double eps, bool verbose) { MITK_INFO(verbose) << "=== mitk::pa::Vector Equal ==="; if (rightHandSide.IsNull() || leftHandSide.IsNull()) { MITK_INFO(verbose) << "Cannot compare nullpointers"; return false; } if (leftHandSide->GetElement(0) - rightHandSide->GetElement(0) > eps) { MITK_INFO(verbose) << "Element[0] not equal"; return false; } if (leftHandSide->GetElement(1) - rightHandSide->GetElement(1) > eps) { MITK_INFO(verbose) << "Element[1] not equal"; return false; } if (leftHandSide->GetElement(2) - rightHandSide->GetElement(2) > eps) { MITK_INFO(verbose) << "Element[2] not equal"; return false; } return true; } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp index ce041d7349..f4f49f6d4b 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp @@ -1,86 +1,98 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPAVesselDrawer.h" #include "mitkPAVesselProperties.h" mitk::pa::VesselDrawer::VesselDrawer() { } mitk::pa::VesselDrawer::~VesselDrawer() { } -void mitk::pa::VesselDrawer::DrawVesselInVolume(Vector::Pointer fromPosition, +void mitk::pa::VesselDrawer::DrawVesselInVolume( VesselProperties::Pointer properties, InSilicoTissueVolume::Pointer volume) { - Vector::Pointer stepSize = properties->GetDirectionVector(); + Vector::Pointer stepDirection = properties->GetDirectionVector(); + Vector::Pointer fromPosition = properties->GetPositionVector()->Clone(); - fromPosition->SetElement(0, fromPosition->GetElement(0) + stepSize->GetElement(0)); - fromPosition->SetElement(1, fromPosition->GetElement(1) + stepSize->GetElement(1)); - fromPosition->SetElement(2, fromPosition->GetElement(2) + stepSize->GetElement(2)); + Vector::Pointer toPosition = properties->GetPositionVector()->Clone(); + Vector::Pointer totalWalkingDistance = stepDirection->Clone(); + totalWalkingDistance->Scale(1.0 / volume->GetSpacing()); - double xPos = fromPosition->GetElement(0); - double yPos = fromPosition->GetElement(1); - double zPos = fromPosition->GetElement(2); + MITK_INFO << "STEP DIRECTION NORM: " << stepDirection->GetNorm(); + MITK_INFO << "TOTAL WALKING DISTANCE NORM: " << totalWalkingDistance->GetNorm(); - if (!volume->IsInsideVolume(xPos, yPos, zPos)) + while (totalWalkingDistance->GetNorm() >= 1) { - 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)) + 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) { - continue; + 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); + } } - 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/test/mitkPhotoacoustic3dVolumeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp index f83ad44e06..bcf813d344 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp @@ -1,224 +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 #include #include "mitkPAVolume.h" #include #include class mitkPhotoacoustic3dVolumeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacoustic3dVolumeTestSuite); MITK_TEST(TestCorrectGetDataAndSetDataBehavior); MITK_TEST(TestCallingConstructorWithNullParameter); MITK_TEST(TestCallingConstructorWithCorrectParameters); MITK_TEST(TestModifyImage); MITK_TEST(TestModifyComplexImage); MITK_TEST(TestConvertToMitkImage); MITK_TEST(TestDeepCopy); MITK_TEST(TestCatchException); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::Volume::Pointer m_Photoacoustic3dVolume; public: void setUp() { } void TestCallingConstructorWithNullParameter() { bool exceptionEncountered = false; try { - m_Photoacoustic3dVolume = mitk::pa::Volume::New(nullptr, 3, 3, 3); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(nullptr, 3, 3, 3, 1); } catch (...) { exceptionEncountered = true; } CPPUNIT_ASSERT(exceptionEncountered); } void TestCallingConstructorWithCorrectParameters() { double* data = new double[1]; data[0] = 3; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1, 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetXDim() == 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetYDim() == 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetZDim() == 1); } void TestModifyImage() { double* data = new double[1]; data[0] = 3; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1, 1); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_Photoacoustic3dVolume->GetData(0, 0, 0)), m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); m_Photoacoustic3dVolume->SetData(17, 0, 0, 0); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 17); } void TestModifyComplexImage() { unsigned int xDim = 4; unsigned int yDim = 7; unsigned int zDim = 12; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = 5; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, xDim, yDim, zDim); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, xDim, yDim, zDim, 1); for (unsigned int z = 0; z < zDim; z++) for (unsigned int y = 0; y < yDim; y++) for (unsigned int x = 0; x < xDim; x++) { CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(x, y, z) == 5); m_Photoacoustic3dVolume->SetData((x + y)*(z + 1), x, y, z); CPPUNIT_ASSERT(abs(m_Photoacoustic3dVolume->GetData(x, y, z) - (x + y)*(z + 1)) < mitk::eps); } } void TestCorrectGetDataAndSetDataBehavior() { unsigned int xDim = 40; unsigned int yDim = 7; unsigned int zDim = 12; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = 0; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, xDim, yDim, zDim); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, xDim, yDim, zDim, 1); for (unsigned int z = 0; z < zDim; z++) for (unsigned int y = 0; y < yDim; y++) for (unsigned int x = 0; x < xDim; x++) { int index = z*xDim*yDim + x*yDim + y; m_Photoacoustic3dVolume->SetData(index, x, y, z); CPPUNIT_ASSERT_MESSAGE(std::to_string(index), m_Photoacoustic3dVolume->GetData(x, y, z) == index); } } void TestConvertToMitkImage() { double* data = new double[6]; data[0] = 3; data[1] = 3; data[2] = 3; data[3] = 3; data[4] = 3; data[5] = 3; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 2, 3); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 2, 3, 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 2) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 0) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 2) == 3); m_Photoacoustic3dVolume->SetData(17, 0, 0, 0); m_Photoacoustic3dVolume->SetData(17, 0, 1, 0); m_Photoacoustic3dVolume->SetData(17, 0, 1, 2); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 17); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 2) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 0) == 17); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 2) == 17); mitk::Image::Pointer mitkImage = m_Photoacoustic3dVolume->AsMitkImage(); CPPUNIT_ASSERT(mitkImage->GetDimensions()[0] == 2); CPPUNIT_ASSERT(mitkImage->GetDimensions()[1] == 1); CPPUNIT_ASSERT(mitkImage->GetDimensions()[2] == 3); mitk::ImageReadAccessor readAccess(mitkImage, mitkImage->GetVolumeData()); double* copyData = (double*)readAccess.GetData(); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[0]), copyData[0] == 17); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[1]), copyData[1] == 17); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[2]), copyData[2] == 3); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[3]), copyData[3] == 3); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[4]), copyData[4] == 3); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[5]), copyData[5] == 17); } void TestDeepCopy() { double* data = new double[1]; data[0] = 3; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1, 1); mitk::pa::Volume::Pointer copiedVolume = m_Photoacoustic3dVolume->DeepCopy(); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetXDim() == copiedVolume->GetXDim()); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetYDim() == copiedVolume->GetYDim()); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetZDim() == copiedVolume->GetZDim()); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); CPPUNIT_ASSERT(copiedVolume->GetData(0, 0, 0) == 3); m_Photoacoustic3dVolume->SetData(17, 0, 0, 0); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 17); CPPUNIT_ASSERT(copiedVolume->GetData(0, 0, 0) == 3); } void AssertIndexException(unsigned int x, unsigned int y, unsigned int z) { bool exceptionCaught = false; try { double thisIsIrrelevant = m_Photoacoustic3dVolume->GetData(x, y, z); thisIsIrrelevant += 1; } catch (...) { exceptionCaught = true; if (exceptionCaught) exceptionCaught = true; } #ifdef _DEBUG CPPUNIT_ASSERT(exceptionCaught); #endif - } + } void TestCatchException() { double* data = new double[1]; data[0] = 3; - m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); + m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1, 1); AssertIndexException(1, 0, 0); AssertIndexException(0, 1, 0); AssertIndexException(0, 0, 1); AssertIndexException(18, 1, 222); } void tearDown() { m_Photoacoustic3dVolume = nullptr; } - }; +}; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacoustic3dVolume) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp index 6d954bba08..5355b82627 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp @@ -1,148 +1,148 @@ /*=================================================================== 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 "mitkPAComposedVolume.h" #include "mitkIOUtil.h" #include "mitkImageReadAccessor.h" #include class mitkPhotoacousticComposedVolumeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticComposedVolumeTestSuite); MITK_TEST(TestCreateAndDestructComposedVolume); MITK_TEST(TestAccessInvalidFluenceComponent); MITK_TEST(TestAccessInvalidFluenceComponentIndex); MITK_TEST(TestAddMultiplePairs); MITK_TEST(TestSortFunctionality); MITK_TEST(TestAccessInvalidFluenceComponentForYOffset); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::ComposedVolume::Pointer m_ComposedVolume; mitk::pa::TissueGeneratorParameters::Pointer m_DefaultParameters; mitk::pa::InSilicoTissueVolume::Pointer m_InSilicoTissueVolume; public: void setUp() { m_DefaultParameters = mitk::pa::TissueGeneratorParameters::New(); m_DefaultParameters->SetXDim(5); m_DefaultParameters->SetYDim(5); m_DefaultParameters->SetZDim(5); m_InSilicoTissueVolume = mitk::pa::InSilicoTissueVolume::New(m_DefaultParameters); m_ComposedVolume = mitk::pa::ComposedVolume::New(m_InSilicoTissueVolume); } mitk::pa::FluenceYOffsetPair::Pointer createFluenceYOffsetPair(double value, double yOffset) { double* data = new double[125]; for (int i = 0; i < 125; ++i) data[i] = value; - mitk::pa::Volume::Pointer volume = mitk::pa::Volume::New(data, 5, 5, 5); + mitk::pa::Volume::Pointer volume = mitk::pa::Volume::New(data, 5, 5, 5, 1); return mitk::pa::FluenceYOffsetPair::New(volume, yOffset); } void TestCreateAndDestructComposedVolume() { CPPUNIT_ASSERT(m_ComposedVolume->GetNumberOfFluenceComponents() == 0); } void TestAccessInvalidFluenceComponent() { bool caughtException = false; try { m_ComposedVolume->GetFluenceValue(0, 0, 0, 0); } catch (mitk::Exception e) { caughtException = true; } CPPUNIT_ASSERT(caughtException); } void TestAddMultiplePairs() { m_ComposedVolume->AddSlice(createFluenceYOffsetPair(0, 0)); CPPUNIT_ASSERT(m_ComposedVolume->GetNumberOfFluenceComponents() == 1); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(1, 1)); CPPUNIT_ASSERT(m_ComposedVolume->GetNumberOfFluenceComponents() == 2); } void TestSortFunctionality() { m_ComposedVolume->AddSlice(createFluenceYOffsetPair(2, 2)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(-1, -1)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(1, 1)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(0, 0)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(-2, -2)); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(0, 0, 2, 0) == 2); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(1, 0, 2, 0) == -1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(2, 0, 2, 0) == 1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(3, 0, 2, 0) == 0); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(4, 0, 2, 0) == -2); m_ComposedVolume->Sort(); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(0, 0, 2, 0) == -2); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(1, 0, 2, 0) == -1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(2, 0, 2, 0) == 0); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(3, 0, 2, 0) == 1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(4, 0, 2, 0) == 2); } void TestAccessInvalidFluenceComponentIndex() { #ifdef _DEBUG m_ComposedVolume->AddSlice(createFluenceYOffsetPair(0, 0)); bool caughtException = false; try { double unusedValue = m_ComposedVolume->GetFluenceValue(0, 1, 2, 300); unusedValue = 0; } catch (mitk::Exception e) { caughtException = true; } CPPUNIT_ASSERT(caughtException); #endif } void TestAccessInvalidFluenceComponentForYOffset() { bool caughtException = false; try { m_ComposedVolume->GetYOffsetForFluenceComponentInPixels(0); } catch (mitk::Exception e) { caughtException = true; } CPPUNIT_ASSERT(caughtException); } void tearDown() { m_ComposedVolume = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticComposedVolume) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp index 37d662697f..5010b20d7d 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp @@ -1,187 +1,187 @@ /*=================================================================== 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 class mitkPhotoacousticIOTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticIOTestSuite); MITK_TEST(testLoadInSilicoTissueNrrdFile); MITK_TEST(testLoad3DVolumeNrrdFile); MITK_TEST(testLoad3DVolumeNrrdFileWithBlur); MITK_TEST(testGetNumberOfNrrdFilesInTestDir); MITK_TEST(testGetChildFoldersFromFolder); MITK_TEST(testLoadFCMs); CPPUNIT_TEST_SUITE_END(); private: const std::string TEST_FOLDER_PATH = "testFiles/"; const std::string TEST_IN_SILICO_VOLUME_PATH = "testInSilicoVolume"; const std::string TEST_3D_Volume_PATH = "test3DVolume"; const std::string TEST_FILE_ENDING = ".nrrd"; const std::string TEST_QUALIFIED_FOLDER_PATH = TEST_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + "/"; const std::string FOLDER_FOLDER = "folder/"; const std::string FCM_PATH = TEST_FOLDER_PATH + "fcms/"; const int NUMBER_OF_NRRD_FILES_IN_TEST_DIR = 2; mitk::pa::TissueGeneratorParameters::Pointer m_VolumeProperties; mitk::pa::InSilicoTissueVolume::Pointer m_TestInSilicoVolume; mitk::pa::Volume::Pointer m_Test3DVolume; public: void setUp() { m_VolumeProperties = createTestVolumeParameters(); m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_VolumeProperties); m_Test3DVolume = createTest3DVolume(5); itk::FileTools::CreateDirectory(TEST_FOLDER_PATH); itk::FileTools::CreateDirectory(TEST_QUALIFIED_FOLDER_PATH); itk::FileTools::CreateDirectory(TEST_FOLDER_PATH + FOLDER_FOLDER + FOLDER_FOLDER); itk::FileTools::CreateDirectory(FCM_PATH); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH)); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_QUALIFIED_FOLDER_PATH)); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH + FOLDER_FOLDER + FOLDER_FOLDER)); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(FCM_PATH)); mitk::IOUtil::Save(m_TestInSilicoVolume->ConvertToMitkImage(), TEST_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + TEST_FILE_ENDING); mitk::IOUtil::Save(m_Test3DVolume->AsMitkImage(), TEST_FOLDER_PATH + TEST_3D_Volume_PATH + TEST_FILE_ENDING); auto yo0 = createTest3DVolume(1)->AsMitkImage(); auto yo1 = createTest3DVolume(2)->AsMitkImage(); yo0->GetPropertyList()->SetStringProperty("y-offset", "0"); yo1->GetPropertyList()->SetStringProperty("y-offset", "1"); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("y-offset")); mitk::IOUtil::Save(yo0, TEST_QUALIFIED_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + "_yo0" + TEST_FILE_ENDING); mitk::IOUtil::Save(yo1, TEST_QUALIFIED_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + "_yo1" + TEST_FILE_ENDING); } mitk::pa::Volume::Pointer createTest3DVolume(double value) { unsigned int xDim = 10; unsigned int yDim = 10; unsigned int zDim = 10; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = value; - return mitk::pa::Volume::New(data, xDim, yDim, zDim); + return mitk::pa::Volume::New(data, xDim, yDim, zDim, 1); } mitk::pa::TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); returnParameters->SetBackgroundAbsorption(0); returnParameters->SetBackgroundScattering(0); returnParameters->SetBackgroundAnisotropy(0); return returnParameters; } void assertEqual(mitk::pa::Volume::Pointer first, mitk::pa::Volume::Pointer second) { CPPUNIT_ASSERT(first->GetXDim() == second->GetXDim()); CPPUNIT_ASSERT(first->GetYDim() == second->GetYDim()); CPPUNIT_ASSERT(first->GetZDim() == second->GetZDim()); for (unsigned int x = 0; x < first->GetXDim(); ++x) for (unsigned int y = 0; y < first->GetYDim(); ++y) for (unsigned int z = 0; z < first->GetZDim(); ++z) { std::string message = "Expected " + std::to_string(first->GetData(x, y, z)) + " but was " + std::to_string(second->GetData(x, y, z)); CPPUNIT_ASSERT_MESSAGE(message, abs(first->GetData(x, y, z) - second->GetData(x, y, z)) < 1e-6); } } void testLoadInSilicoTissueNrrdFile() { auto loadedVolume = mitk::pa::IOUtil::LoadInSilicoTissueVolumeFromNrrdFile(TEST_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + TEST_FILE_ENDING); CPPUNIT_ASSERT(loadedVolume->GetTDim() == m_TestInSilicoVolume->GetTDim()); assertEqual(m_TestInSilicoVolume->GetAbsorptionVolume(), loadedVolume->GetAbsorptionVolume()); assertEqual(m_TestInSilicoVolume->GetScatteringVolume(), loadedVolume->GetScatteringVolume()); assertEqual(m_TestInSilicoVolume->GetAnisotropyVolume(), loadedVolume->GetAnisotropyVolume()); } void testLoad3DVolumeNrrdFile() { auto loadedVolume = mitk::pa::IOUtil::LoadNrrd(TEST_FOLDER_PATH + TEST_3D_Volume_PATH + TEST_FILE_ENDING); assertEqual(loadedVolume, m_Test3DVolume); } void testLoad3DVolumeNrrdFileWithBlur() { auto loadedVolume = mitk::pa::IOUtil::LoadNrrd(TEST_FOLDER_PATH + TEST_3D_Volume_PATH + TEST_FILE_ENDING, 1); assertEqual(loadedVolume, m_Test3DVolume); } void testGetNumberOfNrrdFilesInTestDir() { int numberOfFiles = mitk::pa::IOUtil::GetNumberOfNrrdFilesInDirectory(TEST_FOLDER_PATH); CPPUNIT_ASSERT(numberOfFiles == NUMBER_OF_NRRD_FILES_IN_TEST_DIR); } void testGetChildFoldersFromFolder() { std::vector childFolders = mitk::pa::IOUtil::GetAllChildfoldersFromFolder(TEST_FOLDER_PATH); CPPUNIT_ASSERT(childFolders.size() == 1); CPPUNIT_ASSERT(childFolders[0] == TEST_FOLDER_PATH); childFolders = mitk::pa::IOUtil::GetAllChildfoldersFromFolder(TEST_FOLDER_PATH + FOLDER_FOLDER); MITK_INFO << "ChildFolders: " << childFolders.size(); CPPUNIT_ASSERT(childFolders.size() == 1); CPPUNIT_ASSERT(childFolders[0] == TEST_FOLDER_PATH + FOLDER_FOLDER + "/folder"); } void testLoadFCMs() { auto fcm1 = createTest3DVolume(1); auto fcm2 = createTest3DVolume(2); auto fcm3 = createTest3DVolume(3); auto fcm4 = createTest3DVolume(4); mitk::IOUtil::Save(fcm1->AsMitkImage(), FCM_PATH + "fcm1_p0,0,0FluenceContributionMap.nrrd"); mitk::IOUtil::Save(fcm2->AsMitkImage(), FCM_PATH + "fcm1_p0,0,1FluenceContributionMap.nrrd"); mitk::IOUtil::Save(fcm3->AsMitkImage(), FCM_PATH + "fcm1_p1,0,0FluenceContributionMap.nrrd"); mitk::IOUtil::Save(fcm4->AsMitkImage(), FCM_PATH + "fcm1_p1,0,1FluenceContributionMap.nrrd"); int prog = 0; auto map = mitk::pa::IOUtil::LoadFluenceContributionMaps(FCM_PATH, 0, &prog); assertEqual(fcm1, map[mitk::pa::IOUtil::Position{ 0,0 }]); assertEqual(fcm2, map[mitk::pa::IOUtil::Position{ 0,1 }]); assertEqual(fcm3, map[mitk::pa::IOUtil::Position{ 1,0 }]); assertEqual(fcm4, map[mitk::pa::IOUtil::Position{ 1,1 }]); } void tearDown() { //CPPUNIT_ASSERT_MESSAGE("Resource leak of test files onto hard drive..", itksys::SystemTools::RemoveADirectory(TEST_FOLDER_PATH) == true); } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticIO) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticNoiseGeneratorTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticNoiseGeneratorTest.cpp index 4661aabd0f..0b6a2edddc 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticNoiseGeneratorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticNoiseGeneratorTest.cpp @@ -1,71 +1,71 @@ /*=================================================================== 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 "mitkPAVolume.h" #include "mitkPANoiseGenerator.h" class mitkPhotoacousticNoiseGeneratorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticNoiseGeneratorTestSuite); MITK_TEST(testNoiseGenerator); CPPUNIT_TEST_SUITE_END(); private: public: mitk::pa::Volume::Pointer m_Volume; void setUp() { } void testNoiseGenerator() { int size = 1000 * 100 * 100; double* volume = new double[size]; for (int i = 0; i < size; i++) { volume[i] = 1; } - m_Volume = mitk::pa::Volume::New(volume, 1000, 100, 100); + m_Volume = mitk::pa::Volume::New(volume, 1000, 100, 100, 1); mitk::pa::NoiseGenerator::ApplyNoiseModel(m_Volume, 0.75, 0.1); int negativecounter = 0; for (int i = 0; i < size; i++) { if (m_Volume->GetData()[i] <= 0) { negativecounter++; } } CPPUNIT_ASSERT_EQUAL_MESSAGE("More than one negative: " + std::to_string(negativecounter) + " (" + std::to_string((((double)negativecounter) / size) * 100) + "%)", negativecounter, 0); } void tearDown() { } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticNoiseGenerator) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp index 0306de205f..1f95d64c9b 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp @@ -1,219 +1,240 @@ /*=================================================================== 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 "mitkPAInSilicoTissueVolume.h" #include "mitkPAVector.h" #include "mitkPAVessel.h" #include "mitkIOUtil.h" class mitkPhotoacousticVesselTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVesselTestSuite); MITK_TEST(testEmptyInitializationProperties); MITK_TEST(testWalkInStraightLine); MITK_TEST(testBifurcate); MITK_TEST(testPartialVolume); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::Vessel::Pointer m_TestVessel; mitk::pa::Vessel::CalculateNewVesselPositionCallback m_StraightLine; mitk::pa::Vessel::CalculateNewVesselPositionCallback m_Diverging; mitk::pa::InSilicoTissueVolume::Pointer m_TestInSilicoVolume; mitk::pa::TissueGeneratorParameters::Pointer m_TestVolumeParameters; public: void setUp() { auto params = mitk::pa::VesselProperties::New(); m_TestVessel = mitk::pa::Vessel::New(params); m_StraightLine = &mitk::pa::VesselMeanderStrategy::CalculateNewPositionInStraightLine; m_Diverging = &mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition; m_TestVolumeParameters = createTestVolumeParameters(); m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_TestVolumeParameters); } mitk::pa::TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); returnParameters->SetBackgroundAbsorption(0); returnParameters->SetBackgroundScattering(0); returnParameters->SetBackgroundAnisotropy(0); return returnParameters; } void testEmptyInitializationProperties() { CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == true); } void testPartialVolume() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetDoPartialVolume(true); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(100); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 3)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 3) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 5) - 5.85786438) <= 1e-6); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 3)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 3) - 5.85786438) <= 1e-6); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 3)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 3) - 5.85786438) <= 1e-6); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 5) - 5.85786438) <= 1e-6); } void testWalkInStraightLine() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetDoPartialVolume(false); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(100); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 6, 4)) <= mitk::eps); + mitk::IOUtil::Save(m_TestInSilicoVolume->GetAbsorptionVolume()->AsMitkImage(), "C:/Users/groehl/Desktop/test.nrrd"); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 6)) <= mitk::eps); - - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 6, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 6, 4)) <= mitk::eps); - m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 6)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 6)) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4)) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 6, 4)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5)) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 6)) <= mitk::eps); + m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(3, 4, 4)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 6, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 6, 4)) <= mitk::eps); + + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5)) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 6)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 6)) <= mitk::eps); + + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(3, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(3, 4, 4)) <= mitk::eps); } void testBifurcate() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(1); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == true); std::mt19937 rng; auto bifurcationVessel = m_TestVessel->Bifurcate(&rng); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(bifurcationVessel->CanBifurcate() == false); } void tearDown() { } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVessel) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp index be8b31aa35..36e086dff3 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp @@ -1,347 +1,347 @@ /*=================================================================== 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 "mitkPAInSilicoTissueVolume.h" #include "mitkPATissueGeneratorParameters.h" class mitkPhotoacousticVolumeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVolumeTestSuite); MITK_TEST(TestInitializedTissueContainsOnlyZeros); MITK_TEST(TestConvertedMitkImageContainsOnlyZerosOrAir); MITK_TEST(TestTissueVolumeContainsCorrectAbsorptionNumber); MITK_TEST(TestTissueVolumeContainsCorrectScatteringNumber); MITK_TEST(TestTissueVolumeContainsCorrectAnisotropyNumber); MITK_TEST(testSecondConstructor); MITK_TEST(testCompleteAirVoxelInclusion); MITK_TEST(testHalfAirVoxelInclusion); MITK_TEST(testCompleteAirAndSkinVoxelInclusion); MITK_TEST(testRandomizeCoefficients); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::InSilicoTissueVolume::Pointer m_PhotoacousticVolume; mitk::pa::TissueGeneratorParameters::Pointer m_TissueGeneratorParameters; public: void setUp() { m_TissueGeneratorParameters = mitk::pa::TissueGeneratorParameters::New(); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); } void TestInitializedTissueContainsOnlyZeros() { int dims = 30; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(0); m_TissueGeneratorParameters->SetBackgroundAbsorption(0); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); for (int x = 0; x < dims; x++) { for (int y = 0; y < dims; y++) { for (int z = 0; z < dims; z++) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Every field should be initialized with 0.", abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(x, y, z)) < mitk::eps, true); } } } } void TestConvertedMitkImageContainsOnlyZerosOrAir() { int dims = 30; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); mitk::Image::Pointer testImage = m_PhotoacousticVolume->ConvertToMitkImage(); mitk::ImageReadAccessor imgMemAcc(testImage); double* imagePointer = (double*)imgMemAcc.GetData(); for (int index = 0; index < dims*dims*dims; index++, imagePointer++) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Every voxel in image should be 0.1 or 0.0001.", true, abs(*imagePointer - 0.1) <= mitk::eps || abs(*imagePointer - 0.0001) <= mitk::eps); } } void TestTissueVolumeContainsCorrectAbsorptionNumber() { int dims = 2; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); m_PhotoacousticVolume->SetVolumeValues(0, 0, 0, 0, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 0, 1, 1, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 0, 2, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 1, 3, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 0, 4, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 1, 5, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 0, 6, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 1, 7, 0, 0); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 0.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 1.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 2.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 3.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 1, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 4.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 5.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 6.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 7.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1)); } void TestTissueVolumeContainsCorrectScatteringNumber() { int dims = 2; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); m_PhotoacousticVolume->SetVolumeValues(0, 0, 0, 0, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 0, 1, 0, 1, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 0, 0, 2, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 1, 0, 3, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 0, 0, 4, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 1, 0, 5, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 0, 0, 6, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 1, 0, 7, 0); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 0.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 1.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 2.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 3.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 1, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 4.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 5.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 6.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 7.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1)); } void TestTissueVolumeContainsCorrectAnisotropyNumber() { int dims = 2; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); m_PhotoacousticVolume->SetVolumeValues(0, 0, 0, 0, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 0, 1, 0, 0, 1); m_PhotoacousticVolume->SetVolumeValues(0, 1, 0, 0, 0, 2); m_PhotoacousticVolume->SetVolumeValues(0, 1, 1, 0, 0, 3); m_PhotoacousticVolume->SetVolumeValues(1, 0, 0, 0, 0, 4); m_PhotoacousticVolume->SetVolumeValues(1, 0, 1, 0, 0, 5); m_PhotoacousticVolume->SetVolumeValues(1, 1, 0, 0, 0, 6); m_PhotoacousticVolume->SetVolumeValues(1, 1, 1, 0, 0, 7); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 0.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 1.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 2.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 3.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 1, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 4.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 5.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 6.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 7.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1)); } mitk::pa::Volume::Pointer createTestVolume(double value) { double* data = new double[27]; for (int i = 0; i < 27; ++i) data[i] = value; - return mitk::pa::Volume::New(data, 3, 3, 3); + return mitk::pa::Volume::New(data, 3, 3, 3, 1); } void assertEqual(mitk::pa::Volume::Pointer first, mitk::pa::Volume::Pointer second) { CPPUNIT_ASSERT(first->GetXDim() == second->GetXDim()); CPPUNIT_ASSERT(first->GetYDim() == second->GetYDim()); CPPUNIT_ASSERT(first->GetZDim() == second->GetZDim()); for (unsigned int x = 0; x < first->GetXDim(); ++x) for (unsigned int y = 0; y < first->GetYDim(); ++y) for (unsigned int z = 0; z < first->GetZDim(); ++z) CPPUNIT_ASSERT(abs(first->GetData(x, y, z) - second->GetData(x, y, z)) < mitk::eps); } void testSecondConstructor() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); assertEqual(m_PhotoacousticVolume->GetAbsorptionVolume(), absorption); assertEqual(m_PhotoacousticVolume->GetScatteringVolume(), scattering); assertEqual(m_PhotoacousticVolume->GetAnisotropyVolume(), anisotropy); assertEqual(m_PhotoacousticVolume->GetSegmentationVolume(), segmentation); } void testCompleteAirVoxelInclusion() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(10); m_TissueGeneratorParameters->SetSkinThicknessInMillimeters(0); m_TissueGeneratorParameters->SetAirAbsorption(2); m_TissueGeneratorParameters->SetAirScattering(4); m_TissueGeneratorParameters->SetAirAnisotropy(6); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 6) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 3) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 3) < mitk::eps); } void testRandomizeCoefficients() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(1); mitk::pa::Volume::Pointer anisotropy = createTestVolume(1); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetRandomizePhysicalProperties(true); m_TissueGeneratorParameters->SetRandomizePhysicalPropertiesPercentage(1); m_TissueGeneratorParameters->SetRngSeed(17); m_TissueGeneratorParameters->SetUseRngSeed(true); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 1) < 0.1); } void testCompleteAirAndSkinVoxelInclusion() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(10); m_TissueGeneratorParameters->SetSkinThicknessInMillimeters(10); m_TissueGeneratorParameters->SetAirAbsorption(2); m_TissueGeneratorParameters->SetAirScattering(4); m_TissueGeneratorParameters->SetAirAnisotropy(6); m_TissueGeneratorParameters->SetSkinAbsorption(4); m_TissueGeneratorParameters->SetSkinScattering(8); m_TissueGeneratorParameters->SetSkinAnisotropy(12); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 8) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 6) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 12) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 3) < mitk::eps); } void testHalfAirVoxelInclusion() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(15); m_TissueGeneratorParameters->SetSkinThicknessInMillimeters(0); m_TissueGeneratorParameters->SetAirAbsorption(2); m_TissueGeneratorParameters->SetAirScattering(4); m_TissueGeneratorParameters->SetAirAnisotropy(6); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 1.5) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 3) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 6) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 4.5) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 3) < mitk::eps); } void tearDown() { m_PhotoacousticVolume = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVolume) diff --git a/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp b/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp index 305057568e..b8709dc09d 100644 --- a/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp @@ -1,90 +1,90 @@ ///*=================================================================== //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 class mitkSimulationBatchGeneratorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkSimulationBatchGeneratorTestSuite); MITK_TEST(testGenerateBatchFileString); MITK_TEST(testGenerateBatchFileAndSaveFile); CPPUNIT_TEST_SUITE_END(); private: const std::string TEST_FOLDER_PATH = "testFiles/"; mitk::pa::SimulationBatchGeneratorParameters::Pointer m_Parameters; mitk::pa::Volume::Pointer m_Test3DVolume; public: void setUp() { m_Parameters = mitk::pa::SimulationBatchGeneratorParameters::New(); m_Parameters->SetBinaryPath("binary"); m_Parameters->SetNrrdFilePath(TEST_FOLDER_PATH); m_Parameters->SetNumberOfPhotons(100); m_Parameters->SetTissueName("tissueName"); m_Parameters->SetVolumeIndex(0); m_Parameters->SetYOffsetLowerThresholdInCentimeters(-1); m_Parameters->SetYOffsetUpperThresholdInCentimeters(1); m_Parameters->SetYOffsetStepInCentimeters(0.5); m_Test3DVolume = createTest3DVolume(5); itk::FileTools::CreateDirectory(TEST_FOLDER_PATH); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH)); } mitk::pa::Volume::Pointer createTest3DVolume(double value) { unsigned int xDim = 10; unsigned int yDim = 10; unsigned int zDim = 10; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = value; - return mitk::pa::Volume::New(data, xDim, yDim, zDim); + return mitk::pa::Volume::New(data, xDim, yDim, zDim, 1); } void testGenerateBatchFileString() { std::string batchGenerationString = mitk::pa::SimulationBatchGenerator::CreateBatchSimulationString(m_Parameters); CPPUNIT_ASSERT(!batchGenerationString.empty()); } void testGenerateBatchFileAndSaveFile() { mitk::pa::SimulationBatchGenerator::WriteBatchFileAndSaveTissueVolume(m_Parameters, m_Test3DVolume->AsMitkImage()); CPPUNIT_ASSERT(itksys::SystemTools::FileExists(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000.nrrd")); CPPUNIT_ASSERT(itksys::SystemTools::FileExists(TEST_FOLDER_PATH + "simulate_all.sh") || itksys::SystemTools::FileExists(TEST_FOLDER_PATH + "simulate_all.bat")); CPPUNIT_ASSERT(itksys::SystemTools::FileExists(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000") && itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000")); } void tearDown() { m_Parameters = nullptr; CPPUNIT_ASSERT_MESSAGE("Resource leak of test files onto hard drive..", itksys::SystemTools::RemoveADirectory(TEST_FOLDER_PATH) == true); } }; MITK_TEST_SUITE_REGISTRATION(mitkSimulationBatchGenerator) diff --git a/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp b/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp index 812ce231e8..66e5341de6 100644 --- a/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp @@ -1,187 +1,187 @@ ///*=================================================================== //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 class mitkSlicedVolumeGeneratorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkSlicedVolumeGeneratorTestSuite); MITK_TEST(testConstructorDestructor); MITK_TEST(testGetSlicedFluenceVolume); MITK_TEST(testGetSlicedFluenceVolumeInverse); MITK_TEST(testGetSlicedFluenceVolumeWithPrecorrection); MITK_TEST(testGetSlicedFluenceVolumeWithPrecorrectionInverse); MITK_TEST(testGetSlicedSignalVolume); MITK_TEST(testGetSlicedAbsorptionVolume); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::ComposedVolume::Pointer m_ComposedVolume; mitk::pa::TissueGeneratorParameters::Pointer m_DefaultParameters; mitk::pa::InSilicoTissueVolume::Pointer m_InSilicoTissueVolume; mitk::pa::SlicedVolumeGenerator::Pointer m_SlicedVolumeGenerator; mitk::pa::Volume::Pointer m_PrecorrectionVolume; public: void setUp() { m_SlicedVolumeGenerator = nullptr; m_DefaultParameters = mitk::pa::TissueGeneratorParameters::New(); m_DefaultParameters->SetXDim(3); m_DefaultParameters->SetYDim(3); m_DefaultParameters->SetZDim(3); m_InSilicoTissueVolume = mitk::pa::InSilicoTissueVolume::New(m_DefaultParameters); m_ComposedVolume = mitk::pa::ComposedVolume::New(m_InSilicoTissueVolume); m_ComposedVolume->AddSlice(CreateValidationPair(-1, 1)); m_ComposedVolume->AddSlice(CreateValidationPair(0, 3)); m_ComposedVolume->AddSlice(CreateValidationPair(1, 6)); m_PrecorrectionVolume = CreatePrecorrectionVolume(); } mitk::pa::Volume::Pointer CreatePrecorrectionVolume() { double* data = new double[27]; for (int i = 0; i < 27; ++i) data[i] = 0.5; - return mitk::pa::Volume::New(data, 3, 3, 3); + return mitk::pa::Volume::New(data, 3, 3, 3, 1); } void FillYSliceWith(mitk::pa::Volume::Pointer fluenceVolume, double ySlice, double value) { for (unsigned int x = 0; x < fluenceVolume->GetXDim(); ++x) for (unsigned int z = 0; z < fluenceVolume->GetZDim(); ++z) { fluenceVolume->SetData(value, x, ySlice, z); } } mitk::pa::FluenceYOffsetPair::Pointer CreateValidationPair(double yOffset, int start) { double* data = new double[27]; - mitk::pa::Volume::Pointer fluenceVolume = mitk::pa::Volume::New(data, 3, 3, 3); + mitk::pa::Volume::Pointer fluenceVolume = mitk::pa::Volume::New(data, 3, 3, 3, 1); FillYSliceWith(fluenceVolume, 0, start + 0); FillYSliceWith(fluenceVolume, 1, start + 1); FillYSliceWith(fluenceVolume, 2, start + 2); return mitk::pa::FluenceYOffsetPair::New(fluenceVolume, yOffset); } void AssertYSliceValue(mitk::pa::Volume::Pointer fluenceVolume, double ySlice, double value) { for (unsigned int x = 0; x < fluenceVolume->GetXDim(); ++x) for (unsigned int z = 0; z < fluenceVolume->GetZDim(); ++z) { std::string msg = "Expected: " + std::to_string(value) + " actual: " + std::to_string(fluenceVolume->GetData(x, ySlice, z)); CPPUNIT_ASSERT_MESSAGE(msg, abs(fluenceVolume->GetData(x, ySlice, z) - value) < mitk::eps); } } void testConstructorDestructor() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(0, false, nullptr, false); CPPUNIT_ASSERT(m_SlicedVolumeGenerator.IsNotNull()); } void testGetSlicedFluenceVolume() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1); AssertYSliceValue(slicedFluence, 1, 4); AssertYSliceValue(slicedFluence, 2, 8); } void testGetSlicedFluenceVolumeInverse() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, true); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1); AssertYSliceValue(slicedFluence, 1, 1.0 / 4.0); AssertYSliceValue(slicedFluence, 2, 1.0 / 8.0); } void testGetSlicedFluenceVolumeWithPrecorrection() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, true, m_PrecorrectionVolume, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 2); AssertYSliceValue(slicedFluence, 1, 8); AssertYSliceValue(slicedFluence, 2, 16); } void testGetSlicedFluenceVolumeWithPrecorrectionInverse() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, true, m_PrecorrectionVolume, true); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1.0 / 2); AssertYSliceValue(slicedFluence, 1, 1.0 / 8); AssertYSliceValue(slicedFluence, 2, 1.0 / 16); } void testGetSlicedSignalVolume() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedSignalImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1 * m_DefaultParameters->GetBackgroundAbsorption()); AssertYSliceValue(slicedFluence, 1, 4 * m_DefaultParameters->GetBackgroundAbsorption()); AssertYSliceValue(slicedFluence, 2, 8 * m_DefaultParameters->GetBackgroundAbsorption()); } void testGetSlicedAbsorptionVolume() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedGroundTruthImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, m_DefaultParameters->GetBackgroundAbsorption()); AssertYSliceValue(slicedFluence, 1, m_DefaultParameters->GetBackgroundAbsorption()); AssertYSliceValue(slicedFluence, 2, m_DefaultParameters->GetBackgroundAbsorption()); } void tearDown() { m_SlicedVolumeGenerator = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkSlicedVolumeGenerator)