diff --git a/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h index f82ccf104b..252534f634 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h @@ -1,63 +1,53 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKVESSELDRAWER_H #define MITKVESSELDRAWER_H #include "mitkVector.h" #include "mitkPAVesselMeanderStrategy.h" #include "mitkPAInSilicoTissueVolume.h" #include "mitkPAVector.h" #include "mitkPAVesselProperties.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT VesselDrawer : public itk::LightObject { public: mitkClassMacroItkParent(VesselDrawer, itk::LightObject); itkFactorylessNewMacro(Self); void DrawVesselInVolume(Vector::Pointer toPosition, VesselProperties::Pointer properties, InSilicoTissueVolume::Pointer volume); protected: VesselDrawer(); virtual ~VesselDrawer(); private: }; - /** - * @brief Equal A function comparing two vessel drawer instances for beeing equal - * - * @param rightHandSide A vesseldrawer to be compared - * @param leftHandSide A vesseldrawer to be compared - * @param eps tolarence for comparison. You can use mitk::eps in most cases. - * @param verbose flag indicating if the user wants detailed console output or not. - * @return true, if all subsequent comparisons are true, false otherwise - */ - MITKPHOTOACOUSTICSLIB_EXPORT bool Equal(const VesselDrawer::Pointer leftHandSide, const VesselDrawer::Pointer rightHandSide, double eps, bool verbose); - } + } } #endif // MITKVESSELDRAWER_H diff --git a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp index 2666b3c3c1..5e43ea0e2f 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp @@ -1,352 +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 #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_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; } 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..d94b341b9d 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp @@ -1,137 +1,140 @@ /*=================================================================== 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) { + MITK_INFO << "Initialized by data*"; if (data == nullptr) mitkThrow() << "You may not initialize a mitk::Volume with a nullptr"; m_InternalMitkImage = mitk::Image::New(); unsigned int* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; dimensions[0] = yDim; dimensions[1] = xDim; dimensions[2] = zDim; m_XDim = xDim; m_YDim = yDim; m_ZDim = zDim; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); m_InternalMitkImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); m_InternalMitkImage->SetImportVolume(data, Image::ImportMemoryManagementType::CopyMemory); m_FastAccessDataPointer = GetData(); delete data; } mitk::pa::Volume::Volume(mitk::Image::Pointer image) { + MITK_INFO << "Initialized by mitk::Image"; + if (image.IsNull()) mitkThrow() << "You may not initialize a mitk::Volume with a null reference to an mitk image"; unsigned int* dimensions = image->GetDimensions(); - m_XDim = dimensions[1]; m_YDim = dimensions[0]; + m_XDim = dimensions[1]; m_ZDim = dimensions[2]; m_InternalMitkImage = image; m_FastAccessDataPointer = GetData(); } 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 smartPtr = new mitk::pa::Volume(data, xDim, yDim, zDim); 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()); } 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/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/Utils/mitkPAVesselDrawer.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp index ce041d7349..685e54f193 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp @@ -1,86 +1,73 @@ /*=================================================================== 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, VesselProperties::Pointer properties, InSilicoTissueVolume::Pointer volume) { Vector::Pointer stepSize = properties->GetDirectionVector(); 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)); double xPos = fromPosition->GetElement(0); double yPos = fromPosition->GetElement(1); double zPos = fromPosition->GetElement(2); if (!volume->IsInsideVolume(xPos, yPos, zPos)) { properties->SetRadiusInVoxel(0); return; } double radius = properties->GetRadiusInVoxel(); double ceiledRadius = ceil(radius); for (int x = xPos - ceiledRadius; x <= xPos + ceiledRadius; x += 1) for (int y = yPos - ceiledRadius; y <= yPos + ceiledRadius; y += 1) for (int z = zPos - ceiledRadius; z <= zPos + ceiledRadius; z += 1) { if (!volume->IsInsideVolume(x, y, z)) { continue; } double xDiff = x - xPos; double yDiff = y - yPos; double zDiff = z - zPos; double vectorLengthDiff = radius - sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); if (vectorLengthDiff > 0) { volume->SetVolumeValues(x, y, z, properties->GetAbsorptionCoefficient(), properties->GetScatteringCoefficient(), properties->GetAnisotopyCoefficient(), mitk::pa::InSilicoTissueVolume::SegmentationType::VESSEL); } } } - -bool mitk::pa::Equal(const VesselDrawer::Pointer leftHandSide, const VesselDrawer::Pointer rightHandSide, double eps, bool verbose) -{ - MITK_INFO(verbose) << "=== mitk::pa::Vessel Equal ==="; - - if (rightHandSide.IsNull() || leftHandSide.IsNull()) - { - MITK_INFO(verbose) << "Cannot compare nullpointers"; - return false; - } - - return true; -} diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp index 153a15be90..9a2abb9380 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp @@ -1,220 +1,225 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPAVolumeManipulator.h" #include "mitkPAExceptions.h" #include "mitkPAInSilicoTissueVolume.h" #include "itkResampleImageFilter.h" +#include "itkGaussianInterpolateImageFunction.h" // Includes for image casting between ITK and MITK #include #include +#include + mitk::pa::VolumeManipulator::VolumeManipulator() {} mitk::pa::VolumeManipulator::~VolumeManipulator() {} void mitk::pa::VolumeManipulator::ThresholdImage(Volume::Pointer image, double threshold) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) if (image->GetData(x, y, z) > threshold) image->SetData(1, x, y, z); else image->SetData(0, x, y, z); } void mitk::pa::VolumeManipulator::MultiplyImage(Volume::Pointer image, double factor) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) image->SetData(image->GetData(x, y, z)*factor, x, y, z); } void mitk::pa::VolumeManipulator::Log10Image(Volume::Pointer image) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) { if (image->GetData(x, y, z) < 0) { MITK_ERROR << "Signal was smaller than 0. There is an error in the input image file."; throw InvalidValueException("Signal was smaller than 0. There is an error in the input image file."); } image->SetData(log10(image->GetData(x, y, z)), x, y, z); } } void mitk::pa::VolumeManipulator::RescaleImage(InSilicoTissueVolume::Pointer image, double ratio) { - image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio)); image->SetScatteringVolume(RescaleImage(image->GetScatteringVolume(), ratio)); image->SetAnisotropyVolume(RescaleImage(image->GetAnisotropyVolume(), ratio)); image->SetSegmentationVolume(RescaleImage(image->GetSegmentationVolume(), ratio)); + image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio)); } mitk::pa::Volume::Pointer mitk::pa::VolumeManipulator::RescaleImage(Volume::Pointer image, double ratio) { typedef itk::Image ImageType; typedef itk::ResampleImageFilter FilterType; + typedef itk::GaussianInterpolateImageFunction InterpolatorType; auto input = image->AsMitkImage(); ImageType::Pointer itkInput = ImageType::New(); mitk::CastToItkImage(input, itkInput); ImageType::SizeType outputSize; outputSize[0] = input->GetDimensions()[0] * ratio; outputSize[1] = input->GetDimensions()[1] * ratio; outputSize[2] = input->GetDimensions()[2] * ratio; FilterType::Pointer resampleImageFilter = FilterType::New(); resampleImageFilter->SetInput(itkInput); resampleImageFilter->SetSize(outputSize); - resampleImageFilter->SetOutputSpacing(input->GetGeometry()->GetSpacing()[0] * ratio); + resampleImageFilter->SetInterpolator(InterpolatorType::New()); + resampleImageFilter->SetOutputSpacing(input->GetGeometry()->GetSpacing()[0] / ratio); resampleImageFilter->UpdateLargestPossibleRegion(); ImageType::Pointer output = resampleImageFilter->GetOutput(); mitk::Image::Pointer mitkOutput = mitk::Image::New(); GrabItkImageMemory(output, mitkOutput); return Volume::New(mitkOutput); } /** * @brief Fast 3D Gaussian convolution IIR approximation * @param paVolume * @param sigma * @author Pascal Getreuer * * Copyright (c) 2011, Pascal Getreuer * All rights reserved. * * This program is free software: you can redistribute it and/or modify it * under the terms of the simplified BSD license. * * You should have received a copy of these licenses along with this program. * If not, see . */ void mitk::pa::VolumeManipulator::GaussianBlur3D(mitk::pa::Volume::Pointer paVolume, double sigma) { double* volume = paVolume->GetData(); long width = paVolume->GetYDim(); long height = paVolume->GetXDim(); long depth = paVolume->GetZDim(); const long plane = width*height; const long numel = plane*depth; double lambda, dnu; double nu, boundaryscale, postscale; double *ptr; long i, x, y, z; int step; if (sigma <= 0) return; lambda = (sigma*sigma) / (8.0); dnu = (1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda)) / (2.0*lambda); nu = dnu; boundaryscale = 1.0 / (1.0 - dnu); postscale = pow(dnu / lambda, 12); /* Filter horizontally along each row */ for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { for (step = 0; step < 4; step++) { ptr = volume + width*(y + height*z); ptr[0] *= boundaryscale; /* Filter rightwards */ for (x = 1; x < width; x++) { ptr[x] += nu*ptr[x - 1]; } ptr[x = width - 1] *= boundaryscale; /* Filter leftwards */ for (; x > 0; x--) { ptr[x - 1] += nu*ptr[x]; } } } } /* Filter vertically along each column */ for (z = 0; z < depth; z++) { for (x = 0; x < width; x++) { for (step = 0; step < 4; step++) { ptr = volume + x + plane*z; ptr[0] *= boundaryscale; /* Filter downwards */ for (i = width; i < plane; i += width) { ptr[i] += nu*ptr[i - width]; } ptr[i = plane - width] *= boundaryscale; /* Filter upwards */ for (; i > 0; i -= width) { ptr[i - width] += nu*ptr[i]; } } } } /* Filter along z-dimension */ for (y = 0; y < height; y++) { for (x = 0; x < width; x++) { for (step = 0; step < 4; step++) { ptr = volume + x + width*y; ptr[0] *= boundaryscale; for (i = plane; i < numel; i += plane) { ptr[i] += nu*ptr[i - plane]; } ptr[i = numel - plane] *= boundaryscale; for (; i > 0; i -= plane) { ptr[i - plane] += nu*ptr[i]; } } } } for (i = 0; i < numel; i++) { volume[i] *= postscale; } }