diff --git a/Modules/Core/include/mitkArray.h b/Modules/Core/include/mitkArray.h index 7709718692..c1b86bbd17 100644 --- a/Modules/Core/include/mitkArray.h +++ b/Modules/Core/include/mitkArray.h @@ -1,143 +1,142 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKARRAY_H_ #define MITKARRAY_H_ #include #include "mitkEqual.h" #include "mitkNumericConstants.h" namespace mitk { /** * Methods to copy from itk::FixedArray s (like mitk::Vector and mitk::Point) into ArrayTypes and vice versa. * ArrayTypes here are all types who implement operator[]. * The two templated methods were made free floating so you may specialize them * for your concrete type. */ /** * @brief Copies elements of an array to this Vector * @param[in] array the array whose values will be copied into the Vector. Must be of a type which overrides the [] * operator * @param[out] toArray the FixedArrray (e.g., mitk::Vector or mitk::Point) which should hold the elements of array. * @attention array must be of dimension NVectorDimension! * @attention this method implicitly converts between data types. */ template void FillArray(itk::FixedArray &toArray, const ArrayType &array) { - itk::FixedArray vectorOrPoint; for (unsigned short int var = 0; var < NVectorDimension; ++var) { toArray[var] = array[var]; } } /** * @brief Copies elements of an array to this Vector * @param[in] array the array whose values will be copied into the Vector. Must be of a type which overrides the [] * operator * @return the FixedArray (e.g., mitk::Vector or mitk::Point) which should hold the elements of array. * @attention array must be of dimension NVectorDimension! * @attention this method implicitly converts between data types. */ template itk::FixedArray FillArray(const ArrayType &array) { itk::FixedArray vectorOrPoint; mitk::FillArray(vectorOrPoint, array); return vectorOrPoint; } /** * @brief Copies the elements of this into an array * @param[in] vectorOrPoint the itk::FixedArray which shall be copied into the array. Can e.g. be of type * mitk::Vector * or mitk::Point * @param[out] array the array which will hold the elements. Must be of a type which overrides the [] operator. * @attention array must be of dimension NVectorDimension! * @attention this method implicitly converts between data types. */ template void ToArray(ArrayType &array, const itk::FixedArray &vectorOrPoint) { for (unsigned short int var = 0; var < NVectorDimension; ++var) { array[var] = vectorOrPoint[var]; } } /** * @brief Copies the elements of this into an array * @param[in] vectorOrPoint the itk::FixedArray which shall be copied into the array. Can e.g. be of type * mitk::Vector * or mitk::Point * @return the array which will hold the elements. Must be of a type which overrides the [] operator. * @attention array must be of dimension NVectorDimension! * @attention this method implicitly converts between data types. */ template ArrayType ToArray(const itk::FixedArray &vectorOrPoint) { ArrayType result; mitk::ToArray(result, vectorOrPoint); return result; } // The FillVector3D and FillVector4D methods are implemented for all common array types here template inline void FillVector3D(Tout &out, mitk::ScalarType x, mitk::ScalarType y, mitk::ScalarType z) { out[0] = x; out[1] = y; out[2] = z; } template inline void FillVector4D(Tout &out, mitk::ScalarType x, mitk::ScalarType y, mitk::ScalarType z, mitk::ScalarType t) { out[0] = x; out[1] = y; out[2] = z; out[3] = t; } /** * Compares two ArrayTypes of size size. * ArrayTypes are all Objects/Types who have a [] operator. Pay attention not to set size higher * than the actual size of the ArrayType as this will lead to unexpected results. */ template inline bool EqualArray( TArrayType1 &arrayType1, TArrayType2 &arrayType2, int size, ScalarType eps = mitk::eps, bool verbose = false) { bool isEqual = true; for (int var = 0; var < size; ++var) { isEqual = isEqual && Equal(arrayType1[var], arrayType2[var], eps); } ConditionalOutputOfDifference(arrayType1, arrayType2, eps, verbose, isEqual); return isEqual; } } #endif /* MITKARRAY_H_ */ diff --git a/Modules/Core/test/mitkArbitraryTimeGeometryTest.cpp b/Modules/Core/test/mitkArbitraryTimeGeometryTest.cpp index 4935c25e5d..68f3ed89c5 100644 --- a/Modules/Core/test/mitkArbitraryTimeGeometryTest.cpp +++ b/Modules/Core/test/mitkArbitraryTimeGeometryTest.cpp @@ -1,532 +1,531 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkArbitraryTimeGeometry.h" #include "mitkGeometry3D.h" #include "mitkTestFixture.h" #include "mitkTestingMacros.h" #include class mitkArbitraryTimeGeometryTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkArbitraryTimeGeometryTestSuite); // Test the append method MITK_TEST(CountTimeSteps); MITK_TEST(GetMinimumTimePoint); MITK_TEST(GetMaximumTimePoint); MITK_TEST(GetTimeBounds); MITK_TEST(IsValidTimePoint); MITK_TEST(TimeStepToTimePoint); MITK_TEST(TimePointToTimeStep); MITK_TEST(GetGeometryCloneForTimeStep); MITK_TEST(GetGeometryForTimeStep); MITK_TEST(GetGeometryForTimePoint); MITK_TEST(IsValid); MITK_TEST(Expand); MITK_TEST(ReplaceTimeStepGeometries); MITK_TEST(ClearAllGeometries); MITK_TEST(AppendNewTimeStep); MITK_TEST(HasCollapsedFinalTimeStep); CPPUNIT_TEST_SUITE_END(); private: mitk::Geometry3D::Pointer m_Geometry1; mitk::Geometry3D::Pointer m_Geometry2; mitk::Geometry3D::Pointer m_Geometry3; mitk::Geometry3D::Pointer m_Geometry3_5; mitk::Geometry3D::Pointer m_Geometry4; mitk::Geometry3D::Pointer m_Geometry5; mitk::Geometry3D::Pointer m_InvalidGeometry; mitk::Geometry3D::Pointer m_NewGeometry; mitk::TimePointType m_Geometry1MinTP; mitk::TimePointType m_Geometry2MinTP; mitk::TimePointType m_Geometry3MinTP; mitk::TimePointType m_Geometry3_5MinTP; mitk::TimePointType m_Geometry4MinTP; mitk::TimePointType m_Geometry5MinTP; mitk::TimePointType m_NewGeometryMinTP; mitk::TimePointType m_Geometry1MaxTP; mitk::TimePointType m_Geometry2MaxTP; mitk::TimePointType m_Geometry3MaxTP; mitk::TimePointType m_Geometry3_5MaxTP; mitk::TimePointType m_Geometry4MaxTP; mitk::TimePointType m_Geometry5MaxTP; mitk::TimePointType m_NewGeometryMaxTP; mitk::ArbitraryTimeGeometry::Pointer m_emptyTimeGeometry; mitk::ArbitraryTimeGeometry::Pointer m_initTimeGeometry; mitk::ArbitraryTimeGeometry::Pointer m_12345TimeGeometry; mitk::ArbitraryTimeGeometry::Pointer m_123TimeGeometry; mitk::ArbitraryTimeGeometry::Pointer m_123TimeGeometryWithCollapsedEnd; mitk::ArbitraryTimeGeometry::Pointer m_123TimeGeometryWithCollapsedInterim; public: void setUp() override { - mitk::TimeBounds bounds; m_Geometry1 = mitk::Geometry3D::New(); m_Geometry2 = mitk::Geometry3D::New(); m_Geometry3 = mitk::Geometry3D::New(); m_Geometry3_5 = mitk::Geometry3D::New(); m_Geometry4 = mitk::Geometry3D::New(); m_Geometry5 = mitk::Geometry3D::New(); m_Geometry1MinTP = 1; m_Geometry2MinTP = 2; m_Geometry3MinTP = 3; m_Geometry3_5MinTP = 3.5; m_Geometry4MinTP = 4; m_Geometry5MinTP = 5; m_Geometry1MaxTP = 1.9; m_Geometry2MaxTP = 2.9; m_Geometry3MaxTP = 3.9; m_Geometry3_5MaxTP = 3.9; m_Geometry4MaxTP = 4.9; m_Geometry5MaxTP = 5.9; m_NewGeometry = mitk::Geometry3D::New(); m_NewGeometryMinTP = 20; m_NewGeometryMaxTP = 21.9; mitk::Point3D origin(42); m_NewGeometry->SetOrigin(origin); m_emptyTimeGeometry = mitk::ArbitraryTimeGeometry::New(); m_emptyTimeGeometry->ClearAllGeometries(); m_initTimeGeometry = mitk::ArbitraryTimeGeometry::New(); m_initTimeGeometry->Initialize(); m_12345TimeGeometry = mitk::ArbitraryTimeGeometry::New(); m_12345TimeGeometry->ClearAllGeometries(); m_12345TimeGeometry->AppendNewTimeStep(m_Geometry1, m_Geometry1MinTP, m_Geometry1MaxTP); m_12345TimeGeometry->AppendNewTimeStep(m_Geometry2, m_Geometry2MinTP, m_Geometry2MaxTP); m_12345TimeGeometry->AppendNewTimeStep(m_Geometry3, m_Geometry3MinTP, m_Geometry3MaxTP); m_12345TimeGeometry->AppendNewTimeStep(m_Geometry4, m_Geometry4MinTP, m_Geometry4MaxTP); m_12345TimeGeometry->AppendNewTimeStep(m_Geometry5, m_Geometry5MinTP, m_Geometry5MaxTP); m_123TimeGeometry = mitk::ArbitraryTimeGeometry::New(); m_123TimeGeometry->ClearAllGeometries(); m_123TimeGeometry->AppendNewTimeStep(m_Geometry1, m_Geometry1MinTP, m_Geometry1MaxTP); m_123TimeGeometry->AppendNewTimeStep(m_Geometry2, m_Geometry2MinTP, m_Geometry2MaxTP); m_123TimeGeometry->AppendNewTimeStep(m_Geometry3, m_Geometry3MinTP, m_Geometry3MaxTP); m_123TimeGeometryWithCollapsedEnd = mitk::ArbitraryTimeGeometry::New(); m_123TimeGeometryWithCollapsedEnd->ClearAllGeometries(); m_123TimeGeometryWithCollapsedEnd->AppendNewTimeStep(m_Geometry1, m_Geometry1MinTP, m_Geometry1MaxTP); m_123TimeGeometryWithCollapsedEnd->AppendNewTimeStep(m_Geometry2, m_Geometry2MinTP, m_Geometry2MaxTP); m_123TimeGeometryWithCollapsedEnd->AppendNewTimeStep(m_Geometry3, m_Geometry3MinTP, m_Geometry3MinTP); m_123TimeGeometryWithCollapsedInterim = mitk::ArbitraryTimeGeometry::New(); m_123TimeGeometryWithCollapsedInterim->ClearAllGeometries(); m_123TimeGeometryWithCollapsedInterim->AppendNewTimeStep(m_Geometry1, m_Geometry1MinTP, m_Geometry1MaxTP); m_123TimeGeometryWithCollapsedInterim->AppendNewTimeStep(m_Geometry2, m_Geometry2MinTP, m_Geometry2MinTP); m_123TimeGeometryWithCollapsedInterim->AppendNewTimeStep(m_Geometry3, m_Geometry3MinTP, m_Geometry3MaxTP); } void tearDown() override {} void CountTimeSteps() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->CountTimeSteps() == 0, "Testing CountTimeSteps with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->CountTimeSteps() == 1, "Testing CountTimeSteps with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->CountTimeSteps() == 5, "Testing CountTimeSteps with m_12345TimeGeometry"); } void GetMinimumTimePoint() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMinimumTimePoint() == 0.0, "Testing GetMinimumTimePoint with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetMinimumTimePoint() == 0.0, "Testing GetMinimumTimePoint with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMinimumTimePoint() == 1.0, "Testing GetMinimumTimePoint with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMinimumTimePoint(2) == 0.0, "Testing GetMinimumTimePoint(2) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetMinimumTimePoint(2) == 0.0, "Testing GetMinimumTimePoint(2) with m_initTimeGeometry"); /////////////////////////////////////// // Workarround T27883. See https://phabricator.mitk.org/T27883#219473 for more details. // This workarround should be removed/reevaluated as soon as T28262 is solved and we know // how time geometries should behave in the future! MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMinimumTimePoint(2) == 3.0, "Testing GetMinimumTimePoint(2) with m_12345TimeGeometry"); // Deactivated falling original test // MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMinimumTimePoint(2) == 2.9, // "Testing GetMinimumTimePoint(2) with m_12345TimeGeometry"); // End of workarround for T27883 ////////////////////////////////////// } void GetMaximumTimePoint() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMaximumTimePoint() == 0.0, "Testing GetMaximumTimePoint with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetMaximumTimePoint() == 1.0, "Testing GetMaximumTimePoint with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMaximumTimePoint() == 5.9, "Testing GetMaximumTimePoint with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMaximumTimePoint(2) == 0.0, "Testing GetMaximumTimePoint(2) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetMaximumTimePoint(2) == 0.0, "Testing GetMaximumTimePoint(2) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMaximumTimePoint(2) == 3.9, "Testing GetMaximumTimePoint(2) with m_12345TimeGeometry"); } void GetTimeBounds() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMaximumTimePoint(2) == 0.0, "Testing GetMaximumTimePoint(2) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetMaximumTimePoint(2) == 0.0, "Testing GetMaximumTimePoint(2) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMaximumTimePoint(2) == 3.9, "Testing GetMaximumTimePoint(2) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetTimeBounds()[0] == 0.0, "Testing GetTimeBounds lower part with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetTimeBounds()[0] == 0.0, "Testing GetTimeBounds lower part with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetTimeBounds()[0] == 1.0, "Testing GetTimeBounds lower part with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetTimeBounds()[1] == 0.0, "Testing GetTimeBounds with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetTimeBounds()[1] == 1.0, "Testing GetTimeBounds with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetTimeBounds()[1] == 5.9, "Testing GetTimeBounds with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetTimeBounds(3)[0] == 0.0, "Testing GetTimeBounds(3) lower part with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetTimeBounds(3)[0] == 0.0, "Testing GetTimeBounds(3) lower part with m_initTimeGeometry"); /////////////////////////////////////// // Workarround T27883. See https://phabricator.mitk.org/T27883#219473 for more details. // This workarround should be removed/reevaluated as soon as T28262 is solved and we know // how time geometries should behave in the future! MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetTimeBounds(3)[0] == 4.0, "Testing GetTimeBounds(3) lower part with m_12345TimeGeometry"); // Deactivated falling original test // MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetTimeBounds(3)[0] == 3.9, // "Testing GetTimeBounds(3) lower part with m_12345TimeGeometry"); // End of workarround for T27883 ////////////////////////////////////// MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetTimeBounds(3)[1] == 0.0, "Testing GetTimeBounds(3) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetTimeBounds(3)[1] == 0.0, "Testing GetTimeBounds(3) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetTimeBounds(3)[1] == 4.9, "Testing GetTimeBounds(3) with m_12345TimeGeometry"); } void IsValidTimePoint() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimePoint(-1) == false, "Testing IsValidTimePoint(-1) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimePoint(-1) == false, "Testing IsValidTimePoint(-1) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimePoint(-1) == false, "Testing IsValidTimePoint(-1) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimePoint(0) == false, "Testing IsValidTimePoint(0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimePoint(0) == true, "Testing IsValidTimePoint(0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimePoint(0) == false, "Testing IsValidTimePoint(0) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimePoint(1) == false, "Testing IsValidTimePoint(1) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimePoint(1) == false, "Testing IsValidTimePoint(1) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimePoint(1) == true, "Testing IsValidTimePoint(1) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimePoint(2.5) == false, "Testing IsValidTimePoint(2.5) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimePoint(2.5) == false, "Testing IsValidTimePoint(2.5) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimePoint(2.5) == true, "Testing IsValidTimePoint(2.5) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimePoint(5.89) == false, "Testing IsValidTimePoint(5.89) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimePoint(5.89) == false, "Testing IsValidTimePoint(5.89) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimePoint(5.89) == true, "Testing IsValidTimePoint(5.89) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimePoint(10) == false, "Testing IsValidTimePoint(10) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimePoint(10) == false, "Testing IsValidTimePoint(10) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimePoint(10) == false, "Testing IsValidTimePoint(10) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimeStep(0) == false, "Testing IsValidTimeStep(0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimeStep(0) == true, "Testing IsValidTimeStep(0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimeStep(0) == true, "Testing IsValidTimeStep(0) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimeStep(1) == false, "Testing IsValidTimeStep(1) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimeStep(1) == false, "Testing IsValidTimeStep(1) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimeStep(1) == true, "Testing IsValidTimeStep(1) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValidTimeStep(6) == false, "Testing IsValidTimeStep(6) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValidTimeStep(6) == false, "Testing IsValidTimeStep(6) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValidTimeStep(6) == false, "Testing IsValidTimeStep(6) with m_12345TimeGeometry"); //checked collapsed cases MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometryWithCollapsedInterim->IsValidTimePoint(m_123TimeGeometryWithCollapsedInterim->GetMaximumTimePoint()) == false, "Testing that m_123TimeGeometryWithCollapsedInterim does not inclued the max bound in validity"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometryWithCollapsedEnd->IsValidTimePoint(m_123TimeGeometryWithCollapsedEnd->GetMaximumTimePoint()) == true, "Testing that m_123TimeGeometryWithCollapsedEnd does inclued the max bound in validity, because it has an collapsed final time step. (see also T27259)"); } void TimeStepToTimePoint() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimeStepToTimePoint(0) == 0.0, "Testing TimeStepToTimePoint(0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimeStepToTimePoint(0) == 0.0, "Testing TimeStepToTimePoint(0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimeStepToTimePoint(0) == 1.0, "Testing TimeStepToTimePoint(0) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimeStepToTimePoint(1) == 0.0, "Testing TimeStepToTimePoint(1) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimeStepToTimePoint(1) == 0.0, "Testing TimeStepToTimePoint(1) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimeStepToTimePoint(1) == 2.0, "Testing TimeStepToTimePoint(1) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimeStepToTimePoint(6) == 0.0, "Testing TimeStepToTimePoint(6) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimeStepToTimePoint(6) == 0.0, "Testing TimeStepToTimePoint(6) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimeStepToTimePoint(6) == 0.0, "Testing TimeStepToTimePoint(6) with m_12345TimeGeometry"); } void TimePointToTimeStep() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimePointToTimeStep(0.0) == 0, "Testing TimePointToTimeStep(0.0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimePointToTimeStep(0.0) == 0, "Testing TimePointToTimeStep(0.0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimePointToTimeStep(0.0) == 0, "Testing TimePointToTimeStep(0.0) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimePointToTimeStep(0.5) == 0, "Testing TimePointToTimeStep(0.5) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimePointToTimeStep(0.5) == 0, "Testing TimePointToTimeStep(0.5) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimePointToTimeStep(0.5) == 0, "Testing TimePointToTimeStep(0.5) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimePointToTimeStep(3.5) == 0, "Testing TimePointToTimeStep(3.5) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimePointToTimeStep(3.5) == m_initTimeGeometry->CountTimeSteps(), "Testing TimePointToTimeStep(3.5) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimePointToTimeStep(3.5) == 2, "Testing TimePointToTimeStep(3.5) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->TimePointToTimeStep(5.8) == 0, "Testing TimePointToTimeStep(5.8) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->TimePointToTimeStep(5.8) == m_initTimeGeometry->CountTimeSteps(), "Testing TimePointToTimeStep(5.8) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimePointToTimeStep(5.8) == 4, "Testing TimePointToTimeStep(5.8) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->TimePointToTimeStep(5.9) == m_12345TimeGeometry->CountTimeSteps(), "Testing TimePointToTimeStep(5.9) with m_12345TimeGeometry"); //checked collapsed cases MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometryWithCollapsedInterim->TimePointToTimeStep(m_123TimeGeometryWithCollapsedInterim->GetMaximumTimePoint()) == m_123TimeGeometryWithCollapsedInterim->CountTimeSteps(), "Testing m_123TimeGeometryWithCollapsedInterim does not map the max time poit."); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometryWithCollapsedEnd->TimePointToTimeStep(m_123TimeGeometryWithCollapsedEnd->GetMaximumTimePoint()) == 2, "Testing that m_123TimeGeometryWithCollapsedEnd does map the max bound, because it has an collapsed final time step. (see also T27259)"); } void GetGeometryCloneForTimeStep() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetGeometryCloneForTimeStep(0).IsNull(), "Testing GetGeometryCloneForTimeStep(0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetGeometryCloneForTimeStep(0).IsNotNull(), "Testing GetGeometryCloneForTimeStep(0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetGeometryCloneForTimeStep(0).IsNotNull(), "Testing GetGeometryCloneForTimeStep(0) with m_12345TimeGeometry"); } void GetGeometryForTimeStep() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetGeometryForTimeStep(0).IsNull(), "Testing GetGeometryForTimePoint(0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetGeometryForTimeStep(0).IsNotNull(), "Testing GetGeometryForTimePoint(0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetGeometryForTimeStep(1).IsNull(), "Testing GetGeometryForTimePoint(1) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimeStep(0).GetPointer() == m_Geometry1.GetPointer(), "Testing GetGeometryForTimePoint(0) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimeStep(3).GetPointer() == m_Geometry4.GetPointer(), "Testing GetGeometryForTimePoint(3) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimeStep(4).GetPointer() == m_Geometry5.GetPointer(), "Testing GetGeometryForTimePoint(4) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetGeometryForTimeStep(5).IsNull(), "Testing GetGeometryForTimePoint(5) with m_12345TimeGeometry"); } void GetGeometryForTimePoint() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetGeometryForTimePoint(0).IsNull(), "Testing GetGeometryForTimeStep(0) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetGeometryForTimePoint(0).IsNotNull(), "Testing GetGeometryForTimeStep(0) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetGeometryForTimePoint(0).IsNull(), "Testing GetGeometryForTimeStep(0) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetGeometryForTimePoint(1.5).IsNull(), "Testing GetGeometryForTimeStep(1.5) with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->GetGeometryForTimePoint(1.5).IsNull(), "Testing GetGeometryForTimeStep(1.5) with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimePoint(1.5).GetPointer() == m_Geometry1.GetPointer(), "Testing GetGeometryForTimeStep(1.5) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimePoint(3.5).GetPointer() == m_Geometry3.GetPointer(), "Testing GetGeometryForTimeStep(3.5) with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetGeometryForTimePoint(5.9).IsNull(), "Testing GetGeometryForTimeStep(5.9) with m_12345TimeGeometry"); } void IsValid() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->IsValid() == false, "Testing IsValid() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->IsValid() == true, "Testing IsValid() with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->IsValid() == true, "Testing IsValid() with m_12345TimeGeometry"); } void Expand() { m_12345TimeGeometry->Expand(3); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->CountTimeSteps() == 5, "Testing Expand(3) doesn't change m_12345TimeGeometry"); m_12345TimeGeometry->Expand(7); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->CountTimeSteps() == 7, "Testing Expand(7) with m_12345TimeGeometry"); } void ReplaceTimeStepGeometries() { // Test replace time step geometries m_12345TimeGeometry->ReplaceTimeStepGeometries(m_NewGeometry); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->CountTimeSteps() == 5, "Testing ReplaceTimeStepGeometries() with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimeStep(0)->GetOrigin() == m_NewGeometry->GetOrigin(), "Testing ReplaceTimeStepGeometries(): check if first geometry of m_12345TimeGeometry " "was replaced m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED( m_12345TimeGeometry->GetGeometryForTimeStep(4)->GetOrigin() == m_NewGeometry->GetOrigin(), "Testing ReplaceTimeStepGeometries(): check if last geometry of m_12345TimeGeometry " "was replaced m_12345TimeGeometry"); } void ClearAllGeometries() { // Test clear all geometries m_12345TimeGeometry->ClearAllGeometries(); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->CountTimeSteps() == 0, "Testing ClearAllGeometries() with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMinimumTimePoint() == 0, "Testing ClearAllGeometries() with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->GetMaximumTimePoint() == 0, "Testing ClearAllGeometries() with m_12345TimeGeometry"); } void AppendNewTimeStep() { // Test append MITK_TEST_FOR_EXCEPTION(mitk::Exception, m_123TimeGeometry->AppendNewTimeStep(nullptr, 0, 1)); MITK_TEST_FOR_EXCEPTION(mitk::Exception, m_123TimeGeometry->AppendNewTimeStep(m_Geometry3_5,m_Geometry3_5MinTP,m_Geometry3_5MaxTP)); MITK_TEST_FOR_EXCEPTION(mitk::Exception, m_123TimeGeometry->AppendNewTimeStep(m_Geometry4, m_Geometry4MaxTP, m_Geometry4MinTP)); //valid but inverted bounds m_emptyTimeGeometry->AppendNewTimeStep(m_Geometry4, m_Geometry4MinTP, m_Geometry4MaxTP); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->CountTimeSteps() == 1, "Testing AppendNewTimeStep() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMinimumTimePoint() == 4, "Testing ClearAllGeometries() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->GetMaximumTimePoint() == 4.9, "Testing ClearAllGeometries() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->CountTimeSteps() == 3, "Testing AppendNewTimeStep() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->GetMinimumTimePoint() == 1, "Testing ClearAllGeometries() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->GetMaximumTimePoint() == 3.9, "Testing ClearAllGeometries() with m_emptyTimeGeometry"); m_123TimeGeometry->AppendNewTimeStep(m_Geometry4, m_Geometry4MinTP, m_Geometry4MaxTP); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->CountTimeSteps() == 4, "Testing AppendNewTimeStep() with m_123TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->GetMinimumTimePoint() == 1, "Testing AppendNewTimeStep() with m_123TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->GetMaximumTimePoint() == 4.9, "Testing AppendNewTimeStep() with m_123TimeGeometry"); /////////////////////////////////////// // Workarround T27883. See https://phabricator.mitk.org/T27883#219473 for more details. // This workarround should be removed/reevaluated as soon as T28262 is solved and we know // how time geometries should behave in the future! MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->GetMinimumTimePoint(3) == 4.0, "Testing AppendNewTimeStep() with m_123TimeGeometry"); // Deactivated falling original test // MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometry->GetMinimumTimePoint(3) == 3.9, // "Testing AppendNewTimeStep() with m_123TimeGeometry"); // End of workarround for T27883 ////////////////////////////////////// } void HasCollapsedFinalTimeStep() { MITK_TEST_CONDITION_REQUIRED(m_emptyTimeGeometry->HasCollapsedFinalTimeStep() == false, "Testing HasCollapsedFinalTimeStep() with m_emptyTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_initTimeGeometry->HasCollapsedFinalTimeStep() == false, "Testing HasCollapsedFinalTimeStep() with m_initTimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_12345TimeGeometry->HasCollapsedFinalTimeStep() == false, "Testing HasCollapsedFinalTimeStep() with m_12345TimeGeometry"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometryWithCollapsedEnd->HasCollapsedFinalTimeStep() == true, "Testing HasCollapsedFinalTimeStep() with m_123TimeGeometryWithCollapsedEnd"); MITK_TEST_CONDITION_REQUIRED(m_123TimeGeometryWithCollapsedInterim->HasCollapsedFinalTimeStep() == false, "Testing HasCollapsedFinalTimeStep() with m_123TimeGeometryWithCollapsedInterim"); } }; MITK_TEST_SUITE_REGISTRATION(mitkArbitraryTimeGeometry) diff --git a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx index d301d0d974..706ebb35c1 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx @@ -1,1003 +1,998 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ /*========================================================================= * * Portions of this file are subject to the VTK Toolkit Version 3 copyright. * * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen * * For complete copyright, license and disclaimer of warranty information * please refer to the NOTICE file at the top of the ITK source tree. * *=========================================================================*/ #ifndef __itkMultiGaussianImageSource_hxx #define __itkMultiGaussianImageSource_hxx #include #include #include #include "itkMultiGaussianImageSource.h" #include "itkImageRegionIterator.h" #include "itkObjectFactory.h" #include "itkProgressReporter.h" #include "itkDOMNodeXMLWriter.h" #include namespace itk { /** * */ template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::MultiGaussianImageSource() { //Initial image is 100 wide in each direction. for ( unsigned int i = 0; i < TOutputImage::GetImageDimension(); i++ ) { m_Size[i] = 100; m_Spacing[i] = 1.0; m_Origin[i] = 0.0; m_SphereMidpoint[i] = 0; } m_NumberOfGaussians = 0; m_Radius = 1; m_RadiusStepNumber = 5; m_MeanValue = 0; m_Min = NumericTraits< OutputImagePixelType >::NonpositiveMin(); m_Max = NumericTraits< OutputImagePixelType >::max(); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::~MultiGaussianImageSource() {} //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSize(SizeValueArrayType sizeArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( sizeArray[i] != this->m_Size[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Size[i] = sizeArray[i]; } } } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SizeValueType * MultiGaussianImageSource< TOutputImage > ::GetSize() const { return this->m_Size.GetSize(); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSpacing(SpacingValueArrayType spacingArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( spacingArray[i] != this->m_Spacing[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Spacing[i] = spacingArray[i]; } } } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetOrigin(PointValueArrayType originArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( originArray[i] != this->m_Origin[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Origin[i] = originArray[i]; } } } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::PointValueType * MultiGaussianImageSource< TOutputImage > ::GetOrigin() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_OriginArray[i] = this->m_Origin[i]; } return this->m_OriginArray; } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SpacingValueType * MultiGaussianImageSource< TOutputImage > ::GetSpacing() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_SpacingArray[i] = this->m_Spacing[i]; } return this->m_SpacingArray; } //----------------------------------------------------------------------------------------------------------------------- /** * */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "Max: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Max ) << std::endl; os << indent << "Min: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Min ) << std::endl; os << indent << "Origin: ["; unsigned int ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Origin[ii] << ", "; ++ii; } os << m_Origin[ii] << "]" << std::endl; os << indent << "Spacing: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Spacing[ii] << ", "; ++ii; } os << m_Spacing[ii] << "]" << std::endl; os << indent << "Size: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Size[ii] << ", "; ++ii; } os << m_Size[ii] << "]" << std::endl; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::GetNumberOfGaussians() const { return this->m_NumberOfGaussians; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > typename MultiGaussianImageSource< TOutputImage >::RadiusType MultiGaussianImageSource< TOutputImage > ::GetRadius() const { return this->m_Radius; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRadius( RadiusType radius ) { this->m_Radius = radius; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNumberOfGausssians( unsigned int n ) { this->m_NumberOfGaussians = n; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRegionOfInterest( ItkVectorType roiMin, ItkVectorType roiMax ) { m_RegionOfInterestMax = roiMax; m_RegionOfInterestMin = roiMin; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMaxMeanValue() const { return m_MeanValue; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMaxValueInSphere() const { return m_MaxValueInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetMaxValueIndexInSphere() const { return m_MaxValueIndexInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMinValueInSphere() const { return m_MinValueInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetMinValueIndexInSphere() const { return m_MinValueIndexInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetSphereMidpoint() const { return m_SphereMidpoint; } //----------------------------------------------------------------------------------------------------------------------- /* Calculate and return value of the integral of the gaussian in a cuboid region with the dimension 3: in the x-axis between xMin and xMax and in the y-axis between yMin and yMax and in the z-axis also between zMin and zMax. */ template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtCuboid(double xMin, double xMax, double yMin, double yMax, double zMin, double zMax) { double mean = 0; double summand0, summand1, summand2, value, factor; for(unsigned int n = 0; n < m_NumberOfGaussians; ++n) { summand0 = FunctionPhi((xMax - m_CenterX[n]) / m_SigmaX[n] ) - FunctionPhi((xMin - m_CenterX[n]) / m_SigmaX[n] ); summand1 = FunctionPhi((yMax - m_CenterY[n]) / m_SigmaY[n] ) - FunctionPhi((yMin - m_CenterY[n]) / m_SigmaY[n] ); summand2 = FunctionPhi((zMax - m_CenterZ[n]) / m_SigmaZ[n] ) - FunctionPhi((zMin - m_CenterZ[n]) / m_SigmaZ[n] ); value = summand0 * summand1 * summand2; factor = (m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] ) * pow(2.0 * itk::Math::pi, 1.5 ); mean = mean + factor * value * m_Altitude[n]; } return mean; } //--------------------------------------------------------------------------------------------------------------------- /* Returns the linear interpolation of the values of the standard normal distribution function. This values could be seen in the vector m_NormalDistValues. */ template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::FunctionPhi(double value) { double phiAtValue; //linear interpolation between the values int indexValue = static_cast( 100 * value); if( indexValue > 409 ) { return phiAtValue = 1.0; } else if( indexValue < -409 ) { return phiAtValue = 0.0; } else if( indexValue == 409 ) { return phiAtValue = m_NormalDistValues[409]; } else if( indexValue == -409 ) { return phiAtValue = 1.0 - m_NormalDistValues[409]; } else { bool negative = false; if (indexValue < 0.0) { negative = true; value = -value; } int indexUp = static_cast( 100 * value) + 1; int indexDown = static_cast( 100 * value); double alpha = 100.0 * value - static_cast(indexDown); phiAtValue = (1.0 - alpha) * m_NormalDistValues[indexDown] + alpha * m_NormalDistValues[indexUp] ; if(negative) { phiAtValue = 1.0 - phiAtValue; } return phiAtValue; } } //---------------------------------------------------------------------------------------------------------------------- /* Set the midpoint of the cuboid in a vector m_Midpoints. This cuboids discretise the sphere with the octree method. Set the radius of the cuboid ( = 0.5 * length of the side of the cuboid ) in the vector m_RadiusCuboid. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::InsertPoints( PointType globalCoordinateMidpointCuboid, double cuboidRadius) { typename MapContainerPoints::ElementIdentifier id = m_Midpoints.Size(); m_Midpoints.InsertElement(id, globalCoordinateMidpointCuboid); m_RadiusCuboid.InsertElement(id, cuboidRadius); } //---------------------------------------------------------------------------------------------------------------------- /* This recursive method realise the octree method. It subdivide a cuboid in eight cuboids, when this cuboid crosses the boundary of sphere. If the cuboid is inside the sphere, it calculates the integral and, if uncommented, insert the midpoint of the cuboid in the m_Midpoints vector. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateEdgesInSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, int level ) { double xMin, xMax, yMin, yMax, zMin, zMax; double cuboidRadiusRecursion = cuboidRadius; PointType newMidpoint; int intersect = IntersectTheSphere( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadiusRecursion); if( intersect == 1 ) { if (level < 4) { // cuboid intersect the sphere -> call CalculateEdgesInSphere (this method) for the subdivided cuboid cuboidRadiusRecursion = cuboidRadiusRecursion / 2.0; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { newMidpoint[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadiusRecursion; newMidpoint[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadiusRecursion; newMidpoint[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadiusRecursion; this->CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadiusRecursion, level + 1 ); } } } } // last step of recursion -> on the boundary else { // Calculate the integral and take the half of it (because we are on the boundary) xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // size is in index coordinate -> multiply by the spacing -> global coordinate of the image boundary // yz Plane bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0] * m_Spacing[0]; // xz Plane bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1] * m_Spacing[1]; // xy Plane bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2] * m_Spacing[2]; //check if the boundary of the integral is inside the image if( yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) { //double temp = this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; m_Volume = m_Volume + pow( 2.0 * cuboidRadius, 3.0) / 2.0; } } } else if(intersect == 2) { // cuboid in the sphere // To insert the midpoint and the radius of the cuboid in a vector, so that we can visualise the midpoints, uncomment the next line and the line WriteXMLToTestTheCuboidInsideTheSphere(); // InsertPoints(globalCoordinateMidpointCuboid, cuboidRadius); // Calculate the integral boundary xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // size is in index coordinate -> multiply by the spacing -> global coordinate of the image boundary // yz Plane // bool yzPlaneAtOriginCrossXSection = xMin <= m_Origin[0] && xMax >= m_Origin[0]; // bool yzPlaneAtImageBorderCrossXSection = xMin <= m_Size[0] * m_Spacing[0] && xMax >= m_Size[0] * m_Spacing[0]; bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0] * m_Spacing[0]; // xz Plane // bool xzPlaneAtOriginCrossYSection = yMin <= m_Origin[1] && yMax >= m_Origin[1]; // bool xzPlaneAtImageBorderCrossYSection = yMin <= m_Size[1] * m_Spacing[1] && yMax >= m_Size[1] * m_Spacing[1]; bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1] * m_Spacing[1]; // xy Plane // bool xyPlaneAtOriginCrossZSection = zMin <= m_Origin[2] && zMax >= m_Origin[2]; // bool xyPlaneAtImageBorderCrossZSection = zMin <= m_Size[2] * m_Spacing[2] && zMax >= m_Size[2] * m_Spacing[2]; bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2] * m_Spacing[2]; if( yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) { m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); m_Volume = m_Volume + pow( 2.0 * cuboidRadius, 3.0); } // check if the boundary of the image intersect the cuboid and if yes, change the limits of the cuboid to be only inside the image; therefor we cut the sphere and neglect the part of it outside the image else /* if( // one plane crosses the cuboid ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) || // two plane cross the cuboid (on the image edges possible) ( (yzPlaneAtOriginCrossXSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtImageBorderCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtImageBorderCrossZSection)) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) ) */ { // x-Axis if(xMin <= m_Origin[0] && xMax >= m_Origin[0]) { xMin = m_Origin[0]; }else if(xMin <= m_Size[0] * m_Spacing[0] && xMax >= m_Size[0] * m_Spacing[0]) { xMax = m_Size[0] * m_Spacing[0]; } // y-Axis if(yMin <= m_Origin[1] && yMax >= m_Origin[1]) { yMin = m_Origin[1]; }else if(yMin <= m_Size[1] * m_Spacing[1] && yMax >= m_Size[1] * m_Spacing[1]) { yMax = m_Size[1] * m_Spacing[1]; } // z-Axis if(zMin <= m_Origin[2] && zMax >= m_Origin[2]) { zMin = m_Origin[2]; }else if(zMin <= m_Size[2] * m_Spacing[2] && zMax >= m_Size[2] * m_Spacing[2]) { zMax = m_Size[2] * m_Spacing[2]; } m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin) ; } } } //----------------------------------------------------------------------------------------------------------------------- /* Start the octree recursion in eigth directions for the sphere with midpoint globalCoordinateMidpointSphere and, if uncommented, write this in a file, so that we can visualise it. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateCuboidSegmentationInSphere(PointType globalCoordinateMidpointSphere) { double cuboidRadius = m_Radius * 0.5; PointType newMidpoint; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { newMidpoint[0] = globalCoordinateMidpointSphere[0] + static_cast(i) * cuboidRadius; newMidpoint[1] = globalCoordinateMidpointSphere[1] + static_cast(k) * cuboidRadius; newMidpoint[2] = globalCoordinateMidpointSphere[2] + static_cast(j) * cuboidRadius; CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadius, 0); } } } if(m_WriteMPS) { m_WriteMPS = 0; // uncomment to write an .mps file to visualise the midpoints // std::cout << "Wrote .xml to visualise the midpoints." << std::endl; // WriteXMLToTestTheCuboidInsideTheSphere(); } } //---------------------------------------------------------------------------------------------------------------------- /* This class allows by the method CalculateTheMidpointAndTheMeanValueWithOctree() to find a sphere with a specified radius that has a maximal mean value over all sphere with that radius with midpoint inside or at the boundary of the image. We approximaze the sphere with the octree recursiv method. CalculateTheMidpointAndTheMeanValueWithOctree works as follows: 1. The for-loops traverse the region of interest and assume the current point to be the wanted sphere midpoint. 2. Calculate the mean value for that sphere. 3. Compare with the until-now-found-maximum and take the bigger one. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateTheMidpointAndTheMeanValueWithOctree() { m_MeanValue = 0.0; double meanValueTemp; - PointType midpoint; - //typename MapContainerPoints::ElementIdentifier cuboidNumber = m_Midpoints.Size(); - // SetNormalDistributionValues(); - //double radius; - //double xMin, xMax, yMin, yMax, zMin, zMax; m_WriteMPS = 1; PointType globalCoordinateMidpointSphere; IndexType index; OutputImageRegionType regionOfInterest; IndexType indexR; indexR.SetElement( 0, m_RegionOfInterestMin[0] ); indexR.SetElement( 1, m_RegionOfInterestMin[1] ); indexR.SetElement( 2, m_RegionOfInterestMin[2] ); regionOfInterest.SetIndex(indexR); SizeType sizeROI; sizeROI.SetElement( 0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] + 1); sizeROI.SetElement( 1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] + 1); sizeROI.SetElement( 2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] + 1); regionOfInterest.SetSize(sizeROI); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType regionOfInterestIterator(image, regionOfInterest); for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) { index = regionOfInterestIterator.GetIndex(); image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); m_Volume = 0.0; m_meanValueTemp = 0.0; this->GenerateCuboidSegmentationInSphere(globalCoordinateMidpointSphere); meanValueTemp = m_meanValueTemp / m_Volume; // std::cout << "index: " << index <<" meanvalue: " << meanValueTemp << std::endl; if(meanValueTemp > m_MeanValue) { m_GlobalCoordinate = globalCoordinateMidpointSphere; m_MeanValue = meanValueTemp; m_SphereMidpoint = index; } } } //---------------------------------------------------------------------------------------------------------------------- /* Check if a cuboid intersect the sphere boundary. Returns 2, if the cuboid is inside the sphere; returns 1, if the cuboid intersects the sphere boundary and 0, if the cuboid is out of the sphere. */ template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::IntersectTheSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius) { unsigned int intersect = 1; PointType cuboidEdge; int count = 0; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { cuboidEdge[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadius; cuboidEdge[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadius; cuboidEdge[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadius; if (globalCoordinateMidpointSphere.SquaredEuclideanDistanceTo(cuboidEdge) <= m_Radius * m_Radius) { ++count; } } } } if ( count == 0 ) { // cuboid not in the sphere intersect = 0; } if (count == 8 ) { // cuboid in the sphere intersect = 2; } return intersect; } //---------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtPoint(double x, double y, double z) { //this claculate the mean value in the voxel //integrate over the voxel with midpoint [x, y, z] double summand0, summand1, summand2/*, power*/, value = 0.0, factor; double xMin, xMax, yMin, yMax, zMin, zMax, mean; mean = 0.0; // the for-loop represent the sum of the gaussian function xMin = x - m_Spacing[0] / 2.0; xMax = x + m_Spacing[0] / 2.0; yMin = y - m_Spacing[1] / 2.0; yMax = y + m_Spacing[1] / 2.0; zMin = z - m_Spacing[2] / 2.0; zMax = z + m_Spacing[2] / 2.0; for( unsigned int n = 0; n < m_NumberOfGaussians; ++n ) { summand0 = FunctionPhi( (xMax - m_CenterX[n]) / m_SigmaX[n] ) - FunctionPhi( (xMin - m_CenterX[n]) / m_SigmaX[n] ); summand1 = FunctionPhi( (yMax - m_CenterY[n]) / m_SigmaY[n] ) - FunctionPhi( (yMin - m_CenterY[n]) / m_SigmaY[n] ); summand2 = FunctionPhi( (zMax - m_CenterZ[n]) / m_SigmaZ[n] ) - FunctionPhi( (zMin - m_CenterZ[n]) / m_SigmaZ[n] ); value = summand0 * summand1 * summand2; factor = ( m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] ) * pow( 2.0 * itk::Math::pi, 1.5 ); mean = mean + factor * value * m_Altitude[n]; } value = mean / (m_Spacing[0] * m_Spacing[1] * m_Spacing[2] ); /* //this calculate the value of the gaussian at the midpoint of the voxel: double summand0, summand1, summand2, power, value = 0.0; // the for-loop represent the sum of the gaussian function for(unsigned int n =0; n < m_NumberOfGaussians; ++n) { summand0 = ( x - m_CenterX[n] ) / m_SigmaX[n]; summand1 = ( y - m_CenterY[n] ) / m_SigmaY[n]; summand2 = ( z - m_CenterZ[n] ) / m_SigmaZ[n]; power = summand0 * summand0 + summand1 * summand1 + summand2 * summand2; value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); } */ // std::cout << "X: " << xMin << " " << x << " "<< xMax << " Y: "<< yMin << " " << y << " " << yMax << " Z: "<< zMin << " "<< z << " " << zMax << " value: " << value << std::endl; return value; } //---------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::AddGaussian( VectorType x, VectorType y, VectorType z, VectorType sx, VectorType sy, VectorType sz, VectorType altitude) { for(unsigned int i = 0; i < x.size(); ++i) { m_CenterX.push_back( x[i] ); m_CenterY.push_back( y[i] ); m_CenterZ.push_back( z[i] ); m_SigmaX.push_back( sx[i] ); m_SigmaY.push_back( sy[i] ); m_SigmaZ.push_back( sz[i] ); m_Altitude.push_back(altitude[i]); } } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateOutputInformation() { TOutputImage *output; IndexType index; index.Fill(0); output = this->GetOutput(0); typename TOutputImage::RegionType largestPossibleRegion; largestPossibleRegion.SetSize(this->m_Size); largestPossibleRegion.SetIndex(index); output->SetLargestPossibleRegion(largestPossibleRegion); output->SetSpacing(m_Spacing); output->SetOrigin(m_Origin); } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateData() { itkDebugMacro(<< "Generating a image of scalars "); double valueReal; IndexType index; typename TOutputImage::Pointer image = this->GetOutput(0); image = this->GetOutput(0); image->SetBufferedRegion( image->GetRequestedRegion() ); image->Allocate(); IteratorType imageIt(image, image->GetLargestPossibleRegion()); PointType globalCoordinate; this->SetNormalDistributionValues(); for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) { valueReal = 0.0; index = imageIt.GetIndex(); image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); valueReal = MultiGaussianFunctionValueAtPoint(globalCoordinate[0], globalCoordinate[1] ,globalCoordinate[2]); imageIt.Set(valueReal); } } //----------------------------------------------------------------------------------------------------------------------- /*This method is used to write a .mps file, so that we can visualize the midpoints of the approximated sphere as a scatterplot (for example with MITK Workbench). */ template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::WriteXMLToTestTheCuboidInsideTheSphere() { std::stringstream ss; int numberSummand = 1.0; //write an .mps test file itk::DOMNodeXMLWriter::Pointer xmlWriter; typedef itk::DOMNode::Pointer DOMNodeType; DOMNodeType domXML, domPointSetFile, domFileVersion, domPointSet, domPoint, domId, domX, domY, domZ; xmlWriter = itk::DOMNodeXMLWriter::New(); domXML = itk::DOMNode::New(); domXML->SetName("?xml"); domPointSetFile = itk::DOMNode::New(); domPointSetFile->SetName("point_set_file"); //domFileVersion = itk::DOMNode::New(); //domFileVersion->SetName("file_version"); domPointSet = itk::DOMNode::New(); domPointSet->SetName("point_set"); ss.str(""); ss << 1.0; domXML->SetAttribute("version", ss.str()); domXML->AddChildAtBegin(domPointSetFile); //domPointSetFile -> AddChildAtBegin(domFileVersion); domPointSetFile -> AddChildAtBegin(domPointSet); unsigned int cap = m_Midpoints.Size(); for(unsigned int iter = 0 ; iter < cap; ++iter) { domPoint = itk::DOMNode::New(); domPoint->SetName("point"); domX = itk::DOMNode::New(); domX->SetName("x"); domY = itk::DOMNode::New(); domY->SetName("y"); domZ = itk::DOMNode::New(); domZ->SetName("z"); domId = itk::DOMNode::New(); domId->SetName("id"); ss.str(""); ss << numberSummand; domId->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domId); double scaleFactor = 10.0; PointType point = m_Midpoints.GetElement( numberSummand - 1 ); ss.str(""); ss << point[0] * scaleFactor; domX->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domX); ss.str(""); ss << point[1] * scaleFactor; domY->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domY); ss.str(""); ss << point[2] * scaleFactor; domZ->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domZ); domPointSet -> AddChildAtEnd(domPoint); numberSummand += 1.0; } // .mps (Data) ss.str(""); ss << "C:/temp/CuboidsInTheSphere.mps"; std::string name = ss.str(); char * fileNamePointer = (char*) name.c_str(); xmlWriter->SetFileName( fileNamePointer); xmlWriter->SetInput( domXML ); xmlWriter->Update(); } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateMaxAndMinInSphere() { IndexType index; typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value; m_MaxValueInSphere = std::numeric_limits::min(); m_MinValueInSphere = std::numeric_limits::max(); int radInt, sizeRegion; OutputImageRegionType cuboidRegion; IndexType indexR; SizeType sizeR; int indexRegion, originAsIndex; for( unsigned int i = 0; i < 3; ++i ) { radInt = static_cast(m_Radius/m_Spacing[i]); indexRegion = m_SphereMidpoint[i] - radInt; originAsIndex = static_cast(m_Origin[i]/m_Spacing[i]); if( originAsIndex > indexRegion ) { indexR.SetElement(i, originAsIndex ); } else { indexR.SetElement(i, indexRegion ); } sizeRegion = 2 *radInt + 1; int sizeOutputImage = m_Size[i]; if( (indexR[i] + sizeRegion) > (originAsIndex + sizeOutputImage) ) { std::cout << "Not the entire sphere is in the image!" << std::endl; sizeR.SetElement(i, m_Size[i] - indexRegion ); } else { sizeR.SetElement(i, sizeRegion ); } } cuboidRegion.SetIndex(indexR); cuboidRegion.SetSize(sizeR); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType cuboidRegionOfInterestIterator(image, cuboidRegion); PointType globalCoordinate; for(cuboidRegionOfInterestIterator.GoToBegin(); !cuboidRegionOfInterestIterator.IsAtEnd(); ++cuboidRegionOfInterestIterator) { index = cuboidRegionOfInterestIterator.GetIndex(); if(IsInImage(index)) { image->TransformIndexToPhysicalPoint(cuboidRegionOfInterestIterator.GetIndex(), globalCoordinate); if( m_GlobalCoordinate.EuclideanDistanceTo(globalCoordinate) <= m_Radius ) { value = cuboidRegionOfInterestIterator.Get(); if(m_MaxValueInSphere < value) { m_MaxValueInSphere = value; m_MaxValueIndexInSphere = index; } if(m_MinValueInSphere > value) { m_MinValueInSphere = value; m_MinValueIndexInSphere = index; } } } } } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > bool MultiGaussianImageSource< TOutputImage > ::IsInImage(IndexType index) { bool isInImage = true; int originAsIndex; for( unsigned int i = 0; i < 3; ++i ) { originAsIndex = static_cast(m_Origin[i]/m_Spacing[i]); int sizeOfOutputImage = m_Size[i]; if( index[i] < originAsIndex || index[i] > (originAsIndex + sizeOfOutputImage) ) return false; } return isInImage; } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNormalDistributionValues() { double temp[410] = { 0.5 , 0.50399 , 0.50798, 0.51197, 0.51595, 0.51994, 0.52392, 0.5279, 0.53188, 0.53586, 0.53983, 0.5438, 0.54776, 0.55172, 0.55567, 0.55962, 0.56356, 0.56749, 0.57142, 0.57535, 0.57926, 0.58317, 0.58706, 0.59095, 0.59483 , 0.59871, 0.60257, 0.60642, 0.61026, 0.61409, 0.61791, 0.62172, 0.62552, 0.6293, 0.63307, 0.63683, 0.64058, 0.64431, 0.64803, 0.65173, 0.65542, 0.6591, 0.66276, 0.6664, 0.67003, 0.67364, 0.67724, 0.68082, 0.68439, 0.68793, 0.69146, 0.69497, 0.69847, 0.70194, 0.7054, 0.70884, 0.71226, 0.71566, 0.71904, 0.7224, 0.72575, 0.72907, 0.73237, 0.73565, 0.73891, 0.74215, 0.74537, 0.74857, 0.75175, 0.7549, 0.75804, 0.76115, 0.76424, 0.7673, 0.77035, 0.77337, 0.77637, 0.77935, 0.7823, 0.78524, 0.78814, 0.79103, 0.79389, 0.79673, 0.79955, 0.80234, 0.80511, 0.80785, 0.81057, 0.81327, 0.81594, 0.81859, 0.82121, 0.82381, 0.82639, 0.82894, 0.83147, 0.83398, 0.83646, 0.83891, 0.84134, 0.84375, 0.84614, 0.84849, 0.85083, 0.85314, 0.85543, 0.85769, 0.85993, 0.86214, 0.86433, 0.8665, 0.86864, 0.87076, 0.87286, 0.87493, 0.87698, 0.879, 0.881, 0.88298, 0.88493, 0.88686, 0.88877, 0.89065, 0.89251, 0.89435, 0.89617, 0.89796, 0.89973, 0.90147, 0.9032, 0.9049, 0.90658, 0.90824, 0.90988, 0.91149, 0.91309, 0.91466, 0.91621, 0.91774, 0.91924, 0.92073, 0.9222, 0.92364, 0.92507, 0.92647, 0.92785, 0.92922, 0.93056, 0.93189, 0.93319, 0.93448, 0.93574, 0.93699, 0.93822, 0.93943, 0.94062, 0.94179, 0.94295, 0.94408, 0.9452, 0.9463, 0.94738, 0.94845, 0.9495, 0.95053, 0.95154, 0.95254, 0.95352, 0.95449, 0.95543, 0.95637, 0.95728, 0.95818, 0.95907, 0.95994, 0.9608, 0.96164, 0.96246, 0.96327, 0.96407, 0.96485, 0.96562, 0.96638, 0.96712, 0.96784, 0.96856, 0.96926, 0.96995, 0.97062, 0.97128, 0.97193, 0.97257, 0.9732, 0.97381, 0.97441, 0.975, 0.97558, 0.97615, 0.9767, 0.97725, 0.97778, 0.97831, 0.97882, 0.97932, 0.97982, 0.9803, 0.98077, 0.98124, 0.98169, 0.98214, 0.98257, 0.983, 0.98341, 0.98382, 0.98422, 0.98461, 0.985, 0.98537, 0.98574, 0.9861, 0.98645, 0.98679, 0.98713, 0.98745, 0.98778, 0.98809, 0.9884, 0.9887, 0.98899, 0.98928, 0.98956, 0.98983, 0.9901, 0.99036, 0.99061, 0.99086, 0.99111, 0.99134, 0.99158, 0.9918, 0.99202, 0.99224, 0.99245, 0.99266, 0.99286, 0.99305, 0.99324, 0.99343, 0.99361, 0.99379, 0.99396, 0.99413, 0.9943, 0.99446, 0.99461, 0.99477, 0.99492, 0.99506, 0.9952, 0.99534, 0.99547, 0.9956, 0.99573, 0.99585, 0.99598, 0.99609, 0.99621, 0.99632, 0.99643, 0.99653, 0.99664, 0.99674, 0.99683, 0.99693, 0.99702, 0.99711, 0.9972, 0.99728, 0.99736, 0.99744, 0.99752, 0.9976, 0.99767, 0.99774, 0.99781, 0.99788, 0.99795, 0.99801, 0.99807, 0.99813, 0.99819, 0.99825, 0.99831, 0.99836, 0.99841, 0.99846, 0.99851, 0.99856, 0.99861, 0.99865, 0.99869, 0.99874, 0.99878, 0.99882, 0.99886, 0.99889, 0.99893, 0.99896, 0.999, 0.99903, 0.99906, 0.9991, 0.99913, 0.99916, 0.99918, 0.99921, 0.99924, 0.99926, 0.99929, 0.99931, 0.99934, 0.99936, 0.99938, 0.9994, 0.99942, 0.99944, 0.99946, 0.99948, 0.9995, 0.99952, 0.99953, 0.99955, 0.99957, 0.99958, 0.9996, 0.99961, 0.99962, 0.99964, 0.99965, 0.99966, 0.99968, 0.99969, 0.9997, 0.99971, 0.99972, 0.99973, 0.99974, 0.99975, 0.99976, 0.99977, 0.99978, 0.99978, 0.99979, 0.9998, 0.99981, 0.99981, 0.99982, 0.99983, 0.99983, 0.99984, 0.99985, 0.99985, 0.99986, 0.99986, 0.99987, 0.99987, 0.99988, 0.99988, 0.99989, 0.99989, 0.9999, 0.9999, 0.9999, 0.99991, 0.99991, 0.99992, 0.99992, 0.99992, 0.99992, 0.99993, 0.99993, 0.99993, 0.99994, 0.99994, 0.99994, 0.99994, 0.99995, 0.99995, 0.99995, 0.99995, 0.99995, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99998, 0.99998, 0.99998, 0.99998 }; for(int i=0; i < 410; i++) { m_NormalDistValues[i] = temp[i]; } } } // end namespace itk #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index c35e9d86db..45720f16a3 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,549 +1,548 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void ImageStatisticsCalculator::SetInputImage(const mitk::Image *image) { if (image != m_Image) { m_Image = image; this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator *mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator *mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } mitk::ImageStatisticsContainer* ImageStatisticsCalculator::GetStatistics(LabelIndex label) { if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired(label)) { auto timeGeometry = m_Image->GetTimeGeometry(); // always compute statistics on all timesteps for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); //See T25625: otherwise, the mask is not computed again after setting a different time step m_MaskGenerator->Modified(); m_InternalMask = m_MaskGenerator->GetMask(); if (m_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); imgTimeSel->Update(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } else { // 2) calculate statistics masked AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } } } auto it = m_StatisticContainers.find(label); if (it != m_StatisticContainers.end()) { return (it->second).GetPointer(); } else { mitkThrow() << "unknown label"; return nullptr; } } template void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image *image, const TimeGeometry *timeGeometry, TimeStepType timeStep) { typedef typename itk::Image ImageType; - typedef typename StatisticsImageFilter ImageStatisticsFilterType; + typedef typename mitk::StatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; // reset statistics container if exists ImageStatisticsContainer::Pointer statisticContainerForImage; LabelIndex labelNoMask = 1; auto it = m_StatisticContainers.find(labelNoMask); if (it != m_StatisticContainers.end()) { statisticContainerForImage = it->second; } else { statisticContainerForImage = ImageStatisticsContainer::New(); statisticContainerForImage->SetTimeGeometry(const_cast(timeGeometry)); m_StatisticContainers.emplace(labelNoMask, statisticContainerForImage); } auto statObj = ImageStatisticsContainer::ImageStatisticsObject(); typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: this is single threaded. Implement our own image filter that does this multi threaded // typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = // itk::MinimumMaximumImageCalculator::New(); imgMinMaxFilter->SetImage(image); // imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetInput(image); minMaxFilter->UpdateLargestPossibleRegion(); typename ImageType::PixelType minval = minMaxFilter->GetMin(); typename ImageType::PixelType maxval = minMaxFilter->GetMax(); typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i = 0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i]; maxIndex[i] = tmpMaxIndex[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); // convert m_binSize in m_nBins if necessary unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject &e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } auto voxelVolume = GetVoxelVolume(image); auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels(); auto volume = static_cast(numberOfPixels) * voxelVolume; auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma(); auto rms = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), static_cast(numberOfPixels)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(statisticsFilter->GetMinimum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(statisticsFilter->GetMaximum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); statObj.m_Histogram = statisticsFilter->GetHistogram(); statisticContainerForImage->SetStatisticsForTimeStep(timeStep, statObj); } template double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image *image) const { auto spacing = image->GetSpacing(); double voxelVolume = 1.; for (unsigned int i = 0; i < image->GetImageDimension(); i++) { voxelVolume *= spacing[i]; } return voxelVolume; } template void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename itk::Image *image, const TimeGeometry *timeGeometry, unsigned int timeStep) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef LabelStatisticsImageFilter ImageStatisticsFilterType; typedef MaskUtilities MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; - typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a // 'ignore zuero valued pixels' mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType // unsigned short maskImage = ImageToItkImage(m_InternalMask); } catch (const itk::ExceptionObject &) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this // (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::ConstPointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask // region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue( 1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::unordered_map MapType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::unordered_map nBins; for (LabelPixelType label : relevantLabels) { minVals[label] = static_cast(minMaxFilter->GetMin(label)); maxVals[label] = static_cast(minMaxFilter->GetMax(label)); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } nBins[label] = nBinsForHistogram; } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParameters(nBins, minVals, maxVals); imageStatisticsFilter->Update(); auto labels = imageStatisticsFilter->GetValidLabelValues(); auto it = labels.begin(); while (it != labels.end()) { ImageStatisticsContainer::Pointer statisticContainerForLabelImage; auto labelIt = m_StatisticContainers.find(*it); // reset if statisticContainer already exist if (labelIt != m_StatisticContainers.end()) { statisticContainerForLabelImage = labelIt->second; } // create new statisticContainer else { statisticContainerForLabelImage = ImageStatisticsContainer::New(); statisticContainerForLabelImage->SetTimeGeometry(const_cast(timeGeometry)); // link label (*it) to statisticContainer m_StatisticContainers.emplace(*it, statisticContainerForLabelImage); } ImageStatisticsContainer::ImageStatisticsObject statObj; // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; mitk::Point3D worldCoordinateMin; mitk::Point3D worldCoordinateMax; mitk::Point3D indexCoordinateMin; mitk::Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); minIndex.set_size(3); maxIndex.set_size(3); // for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) for (unsigned int i = 0; i < 3; i++) { minIndex[i] = indexCoordinateMin[i]; maxIndex[i] = indexCoordinateMax[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); auto voxelVolume = GetVoxelVolume(image); auto numberOfVoxels = static_cast(imageStatisticsFilter->GetCount(*it)); auto volume = static_cast(numberOfVoxels) * voxelVolume; auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 auto variance = imageStatisticsFilter->GetSigma(*it) * imageStatisticsFilter->GetSigma(*it); statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(imageStatisticsFilter->GetMinimum(*it))); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(imageStatisticsFilter->GetMaximum(*it))); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it)); statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it); statisticContainerForLabelImage->SetStatisticsForTimeStep(timeStep, statObj); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(LabelIndex label) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); auto it = m_StatisticContainers.find(label); if (it == m_StatisticContainers.end()) { return true; } unsigned long statisticsTimeStamp = it->second->GetMTime(); if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } } // namespace mitk diff --git a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp index 5e8e7529c9..98e8f0dc4e 100644 --- a/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp +++ b/Modules/MatchPointRegistration/include/itkStitchImageFilter.tpp @@ -1,639 +1,638 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef itkStitchImageFilter_hxx #define itkStitchImageFilter_hxx #include "itkStitchImageFilter.h" #include "itkObjectFactory.h" #include "itkIdentityTransform.h" #include "itkProgressReporter.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkImageScanlineIterator.h" #include "itkSpecialCoordinatesImage.h" #include "itkDefaultConvertPixelTraits.h" #include "itkSimpleDataObjectDecorator.h" #include namespace itk { template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::StitchImageFilter() : m_OutputSpacing( 1.0 ), m_OutputOrigin( 0.0 ), m_UseReferenceImage( false ), m_StitchStrategy(StitchStrategy::Mean) { m_Size.Fill( 0 ); m_OutputStartIndex.Fill( 0 ); m_OutputDirection.SetIdentity(); // Pipeline input configuration // implicit input index set: // #1 "ReferenceImage" optional Self::AddOptionalInputName("ReferenceImage"); m_DefaultPixelValue = NumericTraits::ZeroValue( m_DefaultPixelValue ); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(const InputImageType* image) { this->SetInput(0, image, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer(), LinearInterpolatorType::New().GetPointer()); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(unsigned int index, const InputImageType* image) { this->SetInput(index, image, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer(), LinearInterpolatorType::New().GetPointer()); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(unsigned int index, const InputImageType* image, const TransformType* transform) { this->SetInput(index, image, transform, LinearInterpolatorType::New().GetPointer()); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetInput(unsigned int index, const InputImageType* image, const TransformType* transform, InterpolatorType* interpolator) { Superclass::SetInput(index, image); m_Interpolators[image] = interpolator; this->SetTransform(index, transform); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetTransform(unsigned int index, const TransformType* transform) { const auto transformName = this->GetTransformInputName(index); typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; const DecoratorType* oldInput = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(transformName)); if (!oldInput || oldInput->Get() != transform) { typename DecoratorType::Pointer newInput = DecoratorType::New(); // Process object is not const-correct so the const_cast is required here newInput->Set(const_cast(transform)); this->ProcessObject::SetInput(transformName, newInput); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::TransformType* StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetTransform(unsigned int index) const { typedef SimpleDataObjectDecorator< TransformPointerType > DecoratorType; const DecoratorType* input = itkDynamicCastInDebugMode< const DecoratorType* >(this->ProcessObject::GetInput(this->GetTransformInputName(index))); if (nullptr != input) { return input->Get(); } return nullptr; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > const typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >::InterpolatorType* StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetInterpolator(unsigned int index) const { auto input = this->GetInput(index); if (m_Interpolators.find(input) != std::end(m_Interpolators)) { return m_Interpolators[input]; } return nullptr; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputSpacing(const double *spacing) { SpacingType s; for(unsigned int i = 0; i < TOutputImage::ImageDimension; ++i) { s[i] = static_cast< typename SpacingType::ValueType >(spacing[i]); } this->SetOutputSpacing(s); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputOrigin(const double *origin) { OriginPointType p(origin); this->SetOutputOrigin(p); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::SetOutputParametersFromImage(const ImageBaseType *image) { this->SetOutputOrigin ( image->GetOrigin() ); this->SetOutputSpacing ( image->GetSpacing() ); this->SetOutputDirection ( image->GetDirection() ); this->SetOutputStartIndex ( image->GetLargestPossibleRegion().GetIndex() ); this->SetSize ( image->GetLargestPossibleRegion().GetSize() ); } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::BeforeThreadedGenerateData() { this->EnsureInterpolators(); this->EnsureTransforms(); for (const auto& interpolator : m_Interpolators) { interpolator.second->SetInputImage(interpolator.first); } unsigned int nComponents = DefaultConvertPixelTraits::GetNumberOfComponents( m_DefaultPixelValue ); if (nComponents == 0) { PixelComponentType zeroComponent = NumericTraits::ZeroValue(); nComponents = this->GetInput()->GetNumberOfComponentsPerPixel(); NumericTraits::SetLength(m_DefaultPixelValue, nComponents ); for (unsigned int n=0; n void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::AfterThreadedGenerateData() { // Disconnect input image from the interpolator for (auto& interpolator : m_Interpolators) { interpolator.second->SetInputImage(ITK_NULLPTR); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) { if( outputRegionForThread.GetNumberOfPixels() == 0 ) { return; } // Get the output pointers OutputImageType* outputPtr = this->GetOutput(); // Get this input pointers InputImageVectorType inputs = this->GetInputs(); TransformMapType transforms = this->GetTransforms(); std::map lowerIndices; std::map upperIndices; for (const auto& input : inputs) { const auto largestRegion = input->GetLargestPossibleRegion(); lowerIndices[input] = largestRegion.GetIndex(); upperIndices[input] = largestRegion.GetUpperIndex(); } // Create an iterator that will walk the output region for this thread. typedef ImageRegionIteratorWithIndex< OutputImageType > OutputIterator; OutputIterator outIt(outputPtr, outputRegionForThread); // Define a few indices that will be used to translate from an input pixel // to an output pixel PointType outputPoint; // Coordinates of current output pixel PointType inputPoint; // Coordinates of current input pixel ContinuousInputIndexType inputIndex; // Support for progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); // Min/max values of the output pixel type AND these values // represented as the output type of the interpolator const PixelComponentType minValue = NumericTraits< PixelComponentType >::NonpositiveMin(); const PixelComponentType maxValue = NumericTraits< PixelComponentType >::max(); typedef typename InterpolatorType::OutputType OutputType; const ComponentType minOutputValue = static_cast(minValue); const ComponentType maxOutputValue = static_cast(maxValue); // Walk the output region outIt.GoToBegin(); while (!outIt.IsAtEnd()) { // Determine the index of the current output pixel outputPtr->TransformIndexToPhysicalPoint(outIt.GetIndex(), outputPoint); std::vector pixvals; std::vector pixDistance; for (const auto& input : inputs) { // Compute corresponding input pixel position inputPoint = transforms[input]->TransformPoint(outputPoint); const bool isInsideInput = input->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex); // Evaluate input at right position and copy to the output if (m_Interpolators[input]->IsInsideBuffer(inputIndex) && isInsideInput) { OutputType value = m_Interpolators[input]->EvaluateAtContinuousIndex(inputIndex); pixvals.emplace_back(this->CastPixelWithBoundsChecking(value, minOutputValue, maxOutputValue)); - ContinuousInputIndexType indexDistance; const auto spacing = input->GetSpacing(); double minBorderDistance = std::numeric_limits::max(); for (unsigned int i = 0; i < ImageDimension; ++i) { minBorderDistance = std::min(minBorderDistance, std::min(std::abs(lowerIndices[input][i] - inputIndex[i]) * spacing[i], std::abs(upperIndices[input][i] - inputIndex[i]) * spacing[i])); } pixDistance.emplace_back(minBorderDistance); } } if (!pixvals.empty()) { //at least one input provided a value if (StitchStrategy::Mean == m_StitchStrategy) { double sum = std::accumulate(pixvals.begin(), pixvals.end(), 0.0); outIt.Set(sum / pixvals.size()); } else { auto finding = std::max_element(pixDistance.begin(), pixDistance.end()); outIt.Set(pixvals[std::distance(pixDistance.begin(), finding)]); } } else { outIt.Set(m_DefaultPixelValue); // default background value } progress.CompletedPixel(); ++outIt; } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > typename StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::PixelType StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::CastPixelWithBoundsChecking(const InterpolatorOutputType value, const ComponentType minComponent, const ComponentType maxComponent ) const { const unsigned int nComponents = InterpolatorConvertType::GetNumberOfComponents(value); PixelType outputValue; NumericTraits::SetLength( outputValue, nComponents ); for (unsigned int n = 0; n < nComponents; n++) { ComponentType component = InterpolatorConvertType::GetNthComponent( n, value ); if ( component < minComponent ) { PixelConvertType::SetNthComponent( n, outputValue, static_cast( minComponent ) ); } else if ( component > maxComponent ) { PixelConvertType::SetNthComponent( n, outputValue, static_cast( maxComponent ) ); } else { PixelConvertType::SetNthComponent(n, outputValue, static_cast( component ) ); } } return outputValue; } template typename StitchImageFilter::InputImageVectorType StitchImageFilter ::GetInputs() { InputImageVectorType inputs; for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) { auto input = this->GetInput(i); if (nullptr != input) { inputs.push_back(input); } } return inputs; } template typename StitchImageFilter::TransformMapType StitchImageFilter ::GetTransforms() { TransformMapType transforms; for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); ++i) { auto input = this->GetInput(i); auto transform = this->GetTransform(i); transforms[input] = transform; } return transforms; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GenerateInputRequestedRegion() { // Call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); if ( !this->GetInput() ) { return; } // Get pointers to the input auto inputs = this->GetInputs(); for (auto& input : inputs) { InputImagePointer inputPtr = const_cast(input); // Determining the actual input region is non-trivial, especially // when we cannot assume anything about the transform being used. // So we do the easy thing and request the entire input image. // inputPtr->SetRequestedRegionToLargestPossibleRegion(); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GenerateOutputInformation() { // Call the superclass' implementation of this method Superclass::GenerateOutputInformation(); // Get pointers to the input and output OutputImageType *outputPtr = this->GetOutput(); if ( !outputPtr ) { return; } const ReferenceImageBaseType *referenceImage = this->GetReferenceImage(); // Set the size of the output region if ( m_UseReferenceImage && referenceImage ) { outputPtr->SetLargestPossibleRegion( referenceImage->GetLargestPossibleRegion() ); } else { typename TOutputImage::RegionType outputLargestPossibleRegion; outputLargestPossibleRegion.SetSize(m_Size); outputLargestPossibleRegion.SetIndex(m_OutputStartIndex); outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); } // Set spacing and origin if ( m_UseReferenceImage && referenceImage ) { outputPtr->SetSpacing( referenceImage->GetSpacing() ); outputPtr->SetOrigin( referenceImage->GetOrigin() ); outputPtr->SetDirection( referenceImage->GetDirection() ); } else { outputPtr->SetSpacing(m_OutputSpacing); outputPtr->SetOrigin(m_OutputOrigin); outputPtr->SetDirection(m_OutputDirection); } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > ModifiedTimeType StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetMTime(void) const { ModifiedTimeType latestTime = Object::GetMTime(); for (const auto& interpolator : m_Interpolators) { if (interpolator.second.GetPointer()) { if (latestTime < interpolator.second->GetMTime()) { latestTime = interpolator.second->GetMTime(); } } } return latestTime; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "DefaultPixelValue: " << static_cast< typename NumericTraits< PixelType >::PrintType > ( m_DefaultPixelValue ) << std::endl; os << indent << "Size: " << m_Size << std::endl; os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl; os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl; os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl; os << indent << "OutputDirection: " << m_OutputDirection << std::endl; for (const auto& interpolator : m_Interpolators) { os << indent << "Interpolator: " << interpolator.second.GetPointer() << std::endl; } os << indent << "UseReferenceImage: " << ( m_UseReferenceImage ? "On" : "Off" ) << std::endl; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::EnsureTransforms() { const auto inputCount = this->GetNumberOfIndexedInputs(); for (unsigned int i = 0; i < inputCount; ++i) { auto input = this->GetInput(i); if (nullptr == input) { itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); } auto transform = this->GetTransform(i); if (nullptr == transform) { this->SetTransform(i, itk::IdentityTransform< TTransformPrecisionType, ImageDimension>::New().GetPointer()); } } } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > void StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::EnsureInterpolators() { const auto inputCount = this->GetNumberOfIndexedInputs(); InterpolatorMapType newInterpolatorMap; for (unsigned int i = 0; i < inputCount; ++i) { auto input = this->GetInput(i); if (nullptr == input) { itkExceptionMacro(<< "Nth input image is not set (n: " << i << ")."); } if (m_Interpolators[input].IsNull()) { newInterpolatorMap[input] = LinearInterpolatorType::New().GetPointer(); } else { newInterpolatorMap[input] = m_Interpolators[input]; } } m_Interpolators = newInterpolatorMap; } template< typename TInputImage, typename TOutputImage, typename TInterpolatorPrecisionType, typename TTransformPrecisionType > std::string StitchImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType > ::GetTransformInputName(unsigned int index) { return "transform_" + std::to_string(index); } } // end namespace itk #endif