diff --git a/CMakeExternals/ITK.cmake b/CMakeExternals/ITK.cmake index b6ab2759cd..ff402ce770 100644 --- a/CMakeExternals/ITK.cmake +++ b/CMakeExternals/ITK.cmake @@ -1,89 +1,89 @@ #----------------------------------------------------------------------------- # ITK #----------------------------------------------------------------------------- # Sanity checks if(DEFINED ITK_DIR AND NOT EXISTS ${ITK_DIR}) message(FATAL_ERROR "ITK_DIR variable is defined but corresponds to non-existing directory") endif() set(proj ITK) set(proj_DEPENDENCIES GDCM) if(MITK_USE_OpenCV) list(APPEND proj_DEPENDENCIES OpenCV) endif() if(MITK_USE_HDF5) list(APPEND proj_DEPENDENCIES HDF5) endif() set(ITK_DEPENDS ${proj}) if(NOT DEFINED ITK_DIR) set(additional_cmake_args -DUSE_WRAP_ITK:BOOL=OFF) list(APPEND additional_cmake_args -DITKV4_COMPATIBILITY:BOOL=OFF -DITK_LEGACY_REMOVE:BOOL=ON ) if(MITK_USE_OpenCV) list(APPEND additional_cmake_args -DModule_ITKVideoBridgeOpenCV:BOOL=ON -DOpenCV_DIR:PATH=${OpenCV_DIR} "-DCMAKE_CONFIGURATION_TYPES:STRING=Debug$Release" ) endif() # Keep the behaviour of ITK 4.3 which by default turned on ITK Review # see MITK bug #17338 list(APPEND additional_cmake_args -DModule_ITKReview:BOOL=ON -DModule_ITKOpenJPEG:BOOL=ON # for 4.7, the OpenJPEG is needed by review but the variable must be set -DModule_IsotropicWavelets:BOOL=ON ) if(CTEST_USE_LAUNCHERS) list(APPEND additional_cmake_args "-DCMAKE_PROJECT_${proj}_INCLUDE:FILEPATH=${CMAKE_ROOT}/Modules/CTestUseLaunchers.cmake" ) endif() mitk_query_custom_ep_vars() ExternalProject_Add(${proj} LIST_SEPARATOR ${sep} UPDATE_COMMAND "" - GIT_REPOSITORY https://github.com/MITK/ITK.git - GIT_TAG v5.2.1-patched + GIT_REPOSITORY https://github.com/InsightSoftwareConsortium/ITK.git + GIT_TAG v5.3.0 CMAKE_GENERATOR ${gen} CMAKE_GENERATOR_PLATFORM ${gen_platform} CMAKE_ARGS ${ep_common_args} ${additional_cmake_args} -DBUILD_EXAMPLES:BOOL=OFF -DITK_USE_SYSTEM_GDCM:BOOL=ON -DGDCM_DIR:PATH=${GDCM_DIR} -DITK_USE_SYSTEM_HDF5:BOOL=ON -DHDF5_DIR:PATH=${HDF5_DIR} -DModule_GrowCut:BOOL=ON ${${proj}_CUSTOM_CMAKE_ARGS} CMAKE_CACHE_ARGS ${ep_common_cache_args} ${${proj}_CUSTOM_CMAKE_CACHE_ARGS} CMAKE_CACHE_DEFAULT_ARGS ${ep_common_cache_default_args} ${${proj}_CUSTOM_CMAKE_CACHE_DEFAULT_ARGS} DEPENDS ${proj_DEPENDENCIES} ) set(ITK_DIR ${ep_prefix}) mitkFunctionInstallExternalCMakeProject(${proj}) else() mitkMacroEmptyExternalProject(${proj} "${proj_DEPENDENCIES}") endif() diff --git a/CMakeExternals/MatchPoint.cmake b/CMakeExternals/MatchPoint.cmake index cc551e875d..6d988f1658 100644 --- a/CMakeExternals/MatchPoint.cmake +++ b/CMakeExternals/MatchPoint.cmake @@ -1,76 +1,76 @@ #----------------------------------------------------------------------------- # MatchPoint #----------------------------------------------------------------------------- if(MITK_USE_MatchPoint) set(MatchPoint_SOURCE_DIR "" CACHE PATH "Location of the MatchPoint source directory") mark_as_advanced(MatchPoint_SOURCE_DIR) # Sanity checks if(DEFINED MatchPoint_DIR AND NOT EXISTS ${MatchPoint_DIR}) message(FATAL_ERROR "MatchPoint_DIR variable is defined but corresponds to non-existing directory") endif() if(NOT MatchPoint_DIR AND MatchPoint_SOURCE_DIR AND NOT EXISTS ${MatchPoint_SOURCE_DIR}) message(FATAL_ERROR "MatchPoint_SOURCE_DIR variable is defined but corresponds to non-existing directory") endif() set(proj MatchPoint) set(proj_DEPENDENCIES Boost ITK) set(MatchPoint_DEPENDS ${proj}) if(NOT MatchPoint_DIR) set(additional_cmake_args) if(MITK_USE_OpenCV) list(APPEND additional_cmake_args "-DCMAKE_CONFIGURATION_TYPES:STRING=Debug$Release") endif() if(MatchPoint_SOURCE_DIR) set(download_step SOURCE_DIR ${MatchPoint_SOURCE_DIR}) else() set(download_step GIT_REPOSITORY https://github.com/MIC-DKFZ/MatchPoint.git - GIT_TAG e63dfdbb + GIT_TAG a45efdf9 ) endif() string(REPLACE "-DBOOST_ALL_DYN_LINK" "" modified_ep_common_args "${ep_common_args}") ExternalProject_Add(${proj} ${download_step} # INSTALL_COMMAND "${CMAKE_COMMAND} -P cmake_install.cmake" CMAKE_GENERATOR ${gen} CMAKE_GENERATOR_PLATFORM ${gen_platform} CMAKE_ARGS ${modified_ep_common_args} ${additional_cmake_args} -DBUILD_TESTING:BOOL=OFF -DITK_DIR:PATH=${ITK_DIR} #/src/ITK-build "-DBoost_DIR:PATH=${Boost_DIR}" -DMAP_USE_SYSTEM_GDCM:BOOL=ON -DMAP_USE_SYSTEM_HDF5:BOOL=ON -DMAP_DISABLE_ITK_IO_FACTORY_AUTO_REGISTER:BOOL=ON -DMAP_WRAP_Plastimatch:BOOL=ON -DMAP_BUILD_Ontology:BOOL=ON -DMAP_BUILD_Ontology_simple:BOOL=ON -DGDCM_DIR:PATH=${GDCM_DIR} -DHDF5_DIR:PATH=${HDF5_DIR} CMAKE_CACHE_ARGS ${ep_common_cache_args} CMAKE_CACHE_DEFAULT_ARGS ${ep_common_cache_default_args} DEPENDS ${proj_DEPENDENCIES} ) ExternalProject_Get_Property(${proj} binary_dir) set(${proj}_DIR ${binary_dir}) mitkFunctionInstallExternalCMakeProject(${proj}) else() mitkMacroEmptyExternalProject(${proj} "${proj_DEPENDENCIES}") endif() endif(MITK_USE_MatchPoint) diff --git a/Modules/BoundingShape/src/DataManagement/mitkBoundingShapeUtil.cpp b/Modules/BoundingShape/src/DataManagement/mitkBoundingShapeUtil.cpp index f3b18bb3c9..a82b04d219 100644 --- a/Modules/BoundingShape/src/DataManagement/mitkBoundingShapeUtil.cpp +++ b/Modules/BoundingShape/src/DataManagement/mitkBoundingShapeUtil.cpp @@ -1,213 +1,213 @@ /*============================================================================ 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 "mitkBoundingShapeUtil.h" #include "mitkGeometry3D.h" #include "vtkDoubleArray.h" #include "vtkMath.h" #include #include mitk::Handle::Handle() : m_IsActive(false), m_FaceIndices(4), m_Index(0) { m_Position.Fill(0.0); } mitk::Handle::Handle(mitk::Point3D pos, int index, std::vector faceIndices, bool active) : m_IsActive(active), m_Position(pos), m_FaceIndices(faceIndices), m_Index(index) { } mitk::Handle::~Handle() { } bool mitk::Handle::IsActive() { return m_IsActive; } bool mitk::Handle::IsNotActive() { return !m_IsActive; }; void mitk::Handle::SetActive(bool status) { m_IsActive = status; }; void mitk::Handle::SetPosition(mitk::Point3D pos) { m_Position = pos; }; mitk::Point3D mitk::Handle::GetPosition() { return m_Position; }; void mitk::Handle::SetIndex(int index) { m_Index = index; }; int mitk::Handle::GetIndex() { return m_Index; }; std::vector mitk::Handle::GetFaceIndices() { return m_FaceIndices; }; mitk::Point3D mitk::CalcAvgPoint(mitk::Point3D a, mitk::Point3D b) { mitk::Point3D c; c[0] = (a[0] + b[0]) / 2.0; c[1] = (a[1] + b[1]) / 2.0; c[2] = (a[2] + b[2]) / 2.0; return c; } std::vector mitk::GetCornerPoints(mitk::BaseGeometry::Pointer geometry, bool visualizationOffset) { if (geometry == nullptr) mitkThrow() << "Geometry is not valid."; mitk::BoundingBox::ConstPointer boundingBox = geometry->GetBoundingBox(); mitk::Point3D BBmin = boundingBox->GetMinimum(); mitk::Point3D BBmax = boundingBox->GetMaximum(); // use 0.5 offset because the vtkCubeSource is not center pixel based (only for visualization purpose) if (visualizationOffset) { - BBmin -= 0.5; - BBmax -= 0.5; + BBmin -= mitk::Vector(0.5); + BBmax -= mitk::Vector(0.5); } mitk::Point3D p0; p0[0] = BBmin[0]; p0[1] = BBmin[1]; p0[2] = BBmin[2]; // bottom - left - back corner mitk::Point3D p1; p1[0] = BBmin[0]; p1[1] = BBmin[1]; p1[2] = BBmax[2]; // top - left - back corner mitk::Point3D p2; p2[0] = BBmin[0]; p2[1] = BBmax[1]; p2[2] = BBmin[2]; // bottom - left - front corner mitk::Point3D p3; p3[0] = BBmin[0]; p3[1] = BBmax[1]; p3[2] = BBmax[2]; // top - left - front corner mitk::Point3D p4; p4[0] = BBmax[0]; p4[1] = BBmin[1]; p4[2] = BBmin[2]; // bottom - right - back corner mitk::Point3D p5; p5[0] = BBmax[0]; p5[1] = BBmin[1]; p5[2] = BBmax[2]; // top - right - back corner mitk::Point3D p6; p6[0] = BBmax[0]; p6[1] = BBmax[1]; p6[2] = BBmin[2]; // bottom - right - front corner mitk::Point3D p7; p7[0] = BBmax[0]; p7[1] = BBmax[1]; p7[2] = BBmax[2]; // top - right - front corner std::vector cornerPoints; cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p0)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p1)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p2)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p3)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p4)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p5)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p6)); cornerPoints.push_back(geometry->GetIndexToWorldTransform()->TransformPoint(p7)); return cornerPoints; } std::vector mitk::GetHandleIndices(int index) { std::vector faceIndices; faceIndices.resize(4); // +------+ // / /| // +------+ | // | | + // | |/ // +------+ switch (index) { case 0: { faceIndices[0] = 0; faceIndices[1] = 1; faceIndices[2] = 2; faceIndices[3] = 3; } break; case 1: { faceIndices[0] = 4; faceIndices[1] = 5; faceIndices[2] = 6; faceIndices[3] = 7; } break; case 3: { faceIndices[0] = 0; faceIndices[1] = 2; faceIndices[2] = 4; faceIndices[3] = 6; } break; case 2: { faceIndices[0] = 1; faceIndices[1] = 3; faceIndices[2] = 5; faceIndices[3] = 7; } break; case 4: { faceIndices[0] = 0; faceIndices[1] = 1; faceIndices[2] = 4; faceIndices[3] = 5; } break; case 5: { faceIndices[0] = 2; faceIndices[1] = 3; faceIndices[2] = 6; faceIndices[3] = 7; } break; default: { faceIndices[0] = 0; faceIndices[1] = 0; faceIndices[2] = 0; faceIndices[3] = 0; } break; } return faceIndices; } diff --git a/Modules/Classification/CLMiniApps/ManualSegmentationEvaluation.cpp b/Modules/Classification/CLMiniApps/ManualSegmentationEvaluation.cpp index 92ffbb677f..5b6df55277 100644 --- a/Modules/Classification/CLMiniApps/ManualSegmentationEvaluation.cpp +++ b/Modules/Classification/CLMiniApps/ManualSegmentationEvaluation.cpp @@ -1,364 +1,364 @@ /*============================================================================ 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 #include "mitkImage.h" #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include // ITK #include // MITK #include // Classification #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace mitk; std::vector m_FeatureImageVector; void ProcessFeatureImages(const mitk::Image::Pointer & raw_image, const mitk::Image::Pointer & brain_mask) { typedef itk::Image DoubleImageType; typedef itk::Image ShortImageType; typedef itk::ConstNeighborhoodIterator NeighborhoodType; // Neighborhood iterator to access image typedef itk::Functor::NeighborhoodFirstOrderStatistics FunctorType; typedef itk::NeighborhoodFunctorImageFilter FOSFilerType; m_FeatureImageVector.clear(); // RAW m_FeatureImageVector.push_back(raw_image); // GAUSS mitk::Image::Pointer smoothed; mitk::CLUtil::GaussianFilter(raw_image,smoothed,1); m_FeatureImageVector.push_back(smoothed); // Calculate Probability maps (parameters used from literatur) // CSF mitk::Image::Pointer csf_prob = mitk::Image::New(); mitk::CLUtil::ProbabilityMap(smoothed,13.9, 8.3,csf_prob); m_FeatureImageVector.push_back(csf_prob); // Lesion mitk::Image::Pointer les_prob = mitk::Image::New(); mitk::CLUtil::ProbabilityMap(smoothed,59, 11.6,les_prob); m_FeatureImageVector.push_back(les_prob); // Barin (GM/WM) mitk::Image::Pointer brain_prob = mitk::Image::New(); mitk::CLUtil::ProbabilityMap(smoothed,32, 5.6,brain_prob); m_FeatureImageVector.push_back(brain_prob); std::vector FOS_sizes; FOS_sizes.push_back(1); DoubleImageType::Pointer input; ShortImageType::Pointer mask; mitk::CastToItkImage(smoothed, input); mitk::CastToItkImage(brain_mask, mask); for(unsigned int i = 0 ; i < FOS_sizes.size(); i++) { FOSFilerType::Pointer filter = FOSFilerType::New(); filter->SetNeighborhoodSize(FOS_sizes[i]); filter->SetInput(input); filter->SetMask(mask); filter->Update(); FOSFilerType::DataObjectPointerArray array = filter->GetOutputs(); for( unsigned int i = 0; i < FunctorType::OutputCount; i++) { mitk::Image::Pointer featureimage; mitk::CastToMitkImage(dynamic_cast(array[i].GetPointer()),featureimage); m_FeatureImageVector.push_back(featureimage); // AddImageAsDataNode(featureimage,FunctorType::GetFeatureName(i))->SetVisibility(show_nodes); } } { itk::HessianMatrixEigenvalueImageFilter< DoubleImageType >::Pointer filter = itk::HessianMatrixEigenvalueImageFilter< DoubleImageType >::New(); filter->SetInput(input); filter->SetImageMask(mask); filter->SetSigma(3); filter->Update(); mitk::Image::Pointer o1,o2,o3; mitk::CastToMitkImage(filter->GetOutput(0),o1); mitk::CastToMitkImage(filter->GetOutput(1),o2); mitk::CastToMitkImage(filter->GetOutput(2),o3); m_FeatureImageVector.push_back(o1); m_FeatureImageVector.push_back(o2); m_FeatureImageVector.push_back(o3); } { itk::StructureTensorEigenvalueImageFilter< DoubleImageType >::Pointer filter = itk::StructureTensorEigenvalueImageFilter< DoubleImageType >::New(); filter->SetInput(input); filter->SetImageMask(mask); filter->SetInnerScale(1.5); filter->SetOuterScale(3); filter->Update(); mitk::Image::Pointer o1,o2,o3; mitk::CastToMitkImage(filter->GetOutput(0),o1); mitk::CastToMitkImage(filter->GetOutput(1),o2); mitk::CastToMitkImage(filter->GetOutput(2),o3); m_FeatureImageVector.push_back(o1); m_FeatureImageVector.push_back(o2); m_FeatureImageVector.push_back(o3); } { itk::LineHistogramBasedMassImageFilter< DoubleImageType >::Pointer filter = itk::LineHistogramBasedMassImageFilter< DoubleImageType >::New(); filter->SetInput(input); filter->SetImageMask(mask); filter->Update(); mitk::Image::Pointer o1; mitk::CastToMitkImage(filter->GetOutput(0),o1); m_FeatureImageVector.push_back(o1); } } std::vector PointSetToVector(const mitk::PointSet::Pointer & mps) { std::vector result; for(int i = 0 ; i < mps->GetSize(); i++) result.push_back(mps->GetPoint(i)); return result; } int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); // required params parser.addArgument("inputdir", "i", mitkCommandLineParser::Directory, "Input Directory", "Contains input feature files.", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("outputdir", "o", mitkCommandLineParser::Directory, "Output Directory", "Destination of output files.", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("mitkprojectdata", "d", mitkCommandLineParser::File, "original class mask and raw image", "Orig. data.", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("csfmps", "csf", mitkCommandLineParser::File, "CSF Pointset", ".", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("lesmps", "les", mitkCommandLineParser::File, "LES Pointset", ".", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("bramps", "bra", mitkCommandLineParser::File, "BRA Pointset", ".", us::Any(), false, false, false, mitkCommandLineParser::Input); // parser.addArgument("points", "p", mitkCommandLineParser::Int, "Ensure that p points are selected", ".", us::Any(), false); // Miniapp Infos parser.setCategory("Classification Tools"); parser.setTitle("Evaluationtool for Manual-Segmentation"); parser.setDescription("Uses Datacollection to calculate DICE scores for CSF LES BRA"); parser.setContributor("German Cancer Research Center (DKFZ)"); // Params parsing std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inputdir = us::any_cast(parsedArgs["inputdir"]); std::string outputdir = us::any_cast(parsedArgs["outputdir"]); std::string mitkprojectdata = us::any_cast(parsedArgs["mitkprojectdata"]); std::string csf_mps_name = us::any_cast(parsedArgs["csfmps"]); std::string les_mps_name = us::any_cast(parsedArgs["lesmps"]); std::string bra_mps_name = us::any_cast(parsedArgs["bramps"]); mitk::Image::Pointer class_mask_sampled, raw_image, class_mask; mitk::PointSet::Pointer CSF_mps, LES_mps, BRA_mps; // Load from mitk-project auto so = mitk::IOUtil::Load(inputdir + "/" + mitkprojectdata); std::map map; mitk::CLUtil::CountVoxel(dynamic_cast(so[1].GetPointer()), map); raw_image = map.size() <= 7 ? dynamic_cast(so[0].GetPointer()) : dynamic_cast(so[1].GetPointer()); class_mask = map.size() <= 7 ? dynamic_cast(so[1].GetPointer()) : dynamic_cast(so[0].GetPointer()); CSF_mps = mitk::IOUtil::Load(inputdir + "/" + csf_mps_name); LES_mps = mitk::IOUtil::Load(inputdir + "/" + les_mps_name); BRA_mps = mitk::IOUtil::Load(inputdir + "/" + bra_mps_name); unsigned int num_points = CSF_mps->GetSize() + LES_mps->GetSize() + BRA_mps->GetSize(); MITK_INFO << "Found #" << num_points << " points over all classes."; ProcessFeatureImages(raw_image, class_mask); std::map tmpMap; tmpMap[0] = 0; tmpMap[1] = 1; tmpMap[2] = 1; tmpMap[3] = 1; tmpMap[4] = 2; tmpMap[5] = 3; tmpMap[6] = 3; mitk::CLUtil::MergeLabels( class_mask, tmpMap); class_mask_sampled = class_mask->Clone(); itk::Image::Pointer itk_classmask_sampled; mitk::CastToItkImage(class_mask_sampled,itk_classmask_sampled); itk::ImageRegionIteratorWithIndex >::IndexType index; itk::ImageRegionIteratorWithIndex > iit(itk_classmask_sampled,itk_classmask_sampled->GetLargestPossibleRegion()); std::ofstream myfile; myfile.open (inputdir + "/results_3.csv"); Eigen::MatrixXd X_test; unsigned int count_test = 0; mitk::CLUtil::CountVoxel(class_mask, count_test); X_test = Eigen::MatrixXd(count_test, m_FeatureImageVector.size()); unsigned int pos = 0; for( const auto & image : m_FeatureImageVector) { X_test.col(pos) = mitk::CLUtil::Transform(image,class_mask); ++pos; } std::random_device rnd; unsigned int runs = 20; for(unsigned int k = 0 ; k < runs; k++) { auto CSF_vec = PointSetToVector(CSF_mps); auto LES_vec = PointSetToVector(LES_mps); auto BRA_vec = PointSetToVector(BRA_mps); itk_classmask_sampled->FillBuffer(0); // initial draws std::shuffle(CSF_vec.begin(), CSF_vec.end(), std::mt19937(rnd())); class_mask->GetGeometry()->WorldToIndex(CSF_vec.back(),index); iit.SetIndex(index); iit.Set(1); CSF_vec.pop_back(); std::shuffle(LES_vec.begin(), LES_vec.end(), std::mt19937(rnd())); class_mask->GetGeometry()->WorldToIndex(LES_vec.back(),index); iit.SetIndex(index); iit.Set(2); LES_vec.pop_back(); std::shuffle(BRA_vec.begin(), BRA_vec.end(), std::mt19937(rnd())); class_mask->GetGeometry()->WorldToIndex(BRA_vec.back(),index); iit.SetIndex(index); iit.Set(3); BRA_vec.pop_back(); std::stringstream ss; while(!(CSF_vec.empty() && LES_vec.empty() && BRA_vec.empty())) { mitk::CastToMitkImage(itk_classmask_sampled, class_mask_sampled); // Train forest mitk::VigraRandomForestClassifier::Pointer classifier = mitk::VigraRandomForestClassifier::New(); classifier->SetTreeCount(40); classifier->SetSamplesPerTree(0.66); Eigen::MatrixXd X_train; unsigned int count_train = 0; mitk::CLUtil::CountVoxel(class_mask_sampled, count_train); X_train = Eigen::MatrixXd(count_train, m_FeatureImageVector.size() ); unsigned int pos = 0; for( const auto & image : m_FeatureImageVector) { X_train.col(pos) = mitk::CLUtil::Transform(image,class_mask_sampled); ++pos; } Eigen::MatrixXi Y = mitk::CLUtil::Transform(class_mask_sampled,class_mask_sampled); classifier->Train(X_train,Y); Eigen::MatrixXi Y_test = classifier->Predict(X_test); mitk::Image::Pointer result_mask = mitk::CLUtil::Transform(Y_test, class_mask); itk::Image::Pointer itk_result_mask, itk_class_mask; mitk::CastToItkImage(result_mask,itk_result_mask); mitk::CastToItkImage(class_mask, itk_class_mask); itk::LabelOverlapMeasuresImageFilter >::Pointer overlap_filter = itk::LabelOverlapMeasuresImageFilter >::New(); - overlap_filter->SetInput(0,itk_result_mask); - overlap_filter->SetInput(1,itk_class_mask); + overlap_filter->SetSourceImage(itk_result_mask); + overlap_filter->SetTargetImage(itk_class_mask); overlap_filter->Update(); MITK_INFO << "DICE (" << num_points - (CSF_vec.size() + LES_vec.size() + BRA_vec.size()) << "): " << overlap_filter->GetDiceCoefficient(); ss << overlap_filter->GetDiceCoefficient() <<","; // random class selection if(!CSF_vec.empty()) { std::shuffle(CSF_vec.begin(), CSF_vec.end(), std::mt19937(rnd())); class_mask->GetGeometry()->WorldToIndex(CSF_vec.back(),index); iit.SetIndex(index); iit.Set(1); CSF_vec.pop_back(); } if(!LES_vec.empty()) { std::shuffle(LES_vec.begin(), LES_vec.end(), std::mt19937(rnd())); class_mask->GetGeometry()->WorldToIndex(LES_vec.back(),index); iit.SetIndex(index); iit.Set(2); LES_vec.pop_back(); } if(!BRA_vec.empty()) { std::shuffle(BRA_vec.begin(), BRA_vec.end(), std::mt19937(rnd())); class_mask->GetGeometry()->WorldToIndex(BRA_vec.back(),index); iit.SetIndex(index); iit.Set(3); BRA_vec.pop_back(); } } myfile << ss.str() << "\n"; myfile.flush(); } myfile.close(); return EXIT_SUCCESS; } diff --git a/Modules/Core/src/IO/mitkStandardFileLocations.cpp b/Modules/Core/src/IO/mitkStandardFileLocations.cpp index 903f49eefe..12ff940335 100644 --- a/Modules/Core/src/IO/mitkStandardFileLocations.cpp +++ b/Modules/Core/src/IO/mitkStandardFileLocations.cpp @@ -1,239 +1,239 @@ /*============================================================================ 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 #include #include #include #include #include #include mitk::StandardFileLocations::StandardFileLocations() { } mitk::StandardFileLocations::~StandardFileLocations() { } mitk::StandardFileLocations *mitk::StandardFileLocations::GetInstance() { static StandardFileLocations::Pointer m_Instance = nullptr; if (m_Instance.IsNull()) m_Instance = StandardFileLocations::New(); return m_Instance; } void mitk::StandardFileLocations::AddDirectoryForSearch(const char *dir, bool insertInFrontOfSearchList) { // Do nothing if directory is already included into search list (TODO more clever: search only once!) FileSearchVectorType::iterator iter; if (m_SearchDirectories.size() > 0) { iter = std::find(m_SearchDirectories.begin(), m_SearchDirectories.end(), std::string(dir)); if (iter != m_SearchDirectories.end()) return; } // insert dir into queue if (insertInFrontOfSearchList) { auto it = m_SearchDirectories.begin(); m_SearchDirectories.insert(it, std::string(dir)); } else m_SearchDirectories.push_back(std::string(dir)); } void mitk::StandardFileLocations::RemoveDirectoryForSearch(const char *dir) { FileSearchVectorType::iterator it; // background layers if (m_SearchDirectories.size() > 0) { it = std::find(m_SearchDirectories.begin(), m_SearchDirectories.end(), std::string(dir)); if (it != m_SearchDirectories.end()) { m_SearchDirectories.erase(it); return; } } } std::string mitk::StandardFileLocations::SearchDirectoriesForFile(const char *filename) { FileSearchVectorType::iterator it; for (it = m_SearchDirectories.begin(); it != m_SearchDirectories.end(); ++it) { std::string currDir = (*it); // Perhaps append "/" before appending filename if (currDir.find_last_of("\\") + 1 != currDir.size() || currDir.find_last_of("/") + 1 != currDir.size()) currDir += "/"; // Append filename currDir += filename; // Perhaps remove "/" after filename if (currDir.find_last_of("\\") + 1 == currDir.size() || currDir.find_last_of("/") + 1 == currDir.size()) currDir.erase(currDir.size() - 1, currDir.size()); // convert to OS dependent path schema currDir = Utf8Util::Utf8ToLocal8Bit(itksys::SystemTools::ConvertToOutputPath(Utf8Util::Local8BitToUtf8(currDir).c_str())); // On windows systems, the ConvertToOutputPath method quotes paths that contain empty spaces. // These quotes are not expected by the FileExists method and therefore removed, if existing. if (currDir.find_last_of("\"") + 1 == currDir.size()) currDir.erase(currDir.size() - 1, currDir.size()); if (currDir.find_last_of("\"") == 0) currDir.erase(0, 1); // Return first found path if (itksys::SystemTools::FileExists(Utf8Util::Local8BitToUtf8(currDir).c_str())) return currDir; } return std::string(""); } std::string mitk::StandardFileLocations::FindFile(const char *filename, const char *pathInSourceDir) { std::string directoryPath; // 1. look for MITKCONF environment variable const char *mitkConf = itksys::SystemTools::GetEnv("MITKCONF"); if (mitkConf != nullptr) AddDirectoryForSearch(mitkConf, false); // 2. use .mitk-subdirectory in home directory of the user #if defined(_WIN32) && !defined(__CYGWIN__) const char *homeDrive = itksys::SystemTools::GetEnv("HOMEDRIVE"); const char *homePath = itksys::SystemTools::GetEnv("HOMEPATH"); if ((homeDrive != nullptr) || (homePath != nullptr)) { directoryPath = homeDrive; directoryPath += homePath; directoryPath += "/.mitk/"; AddDirectoryForSearch(directoryPath.c_str(), false); } #else const char *homeDirectory = itksys::SystemTools::GetEnv("HOME"); if (homeDirectory != nullptr) { directoryPath = homeDirectory; directoryPath += "/.mitk/"; AddDirectoryForSearch(directoryPath.c_str(), false); } #endif // defined(_WIN32) && !defined(__CYGWIN__) // 3. look in the current working directory directoryPath = ""; AddDirectoryForSearch(directoryPath.c_str()); directoryPath = Utf8Util::Utf8ToLocal8Bit(itksys::SystemTools::GetCurrentWorkingDirectory()); AddDirectoryForSearch(directoryPath.c_str(), false); std::string directoryBinPath = directoryPath + "/bin"; AddDirectoryForSearch(directoryBinPath.c_str(), false); // 4. use a source tree location from compile time directoryPath = MITK_ROOT; if (pathInSourceDir) { directoryPath += pathInSourceDir; } directoryPath += '/'; AddDirectoryForSearch(directoryPath.c_str(), false); return SearchDirectoriesForFile(filename); } std::string mitk::StandardFileLocations::GetOptionDirectory() { const char *mitkoptions = itksys::SystemTools::GetEnv("MITKOPTIONS"); std::string optionsDirectory; if (mitkoptions != nullptr) { // 1. look for MITKOPTIONS environment variable optionsDirectory = mitkoptions; optionsDirectory += "/"; } else { // 2. use .mitk-subdirectory in home directory of the user std::string homeDirectory; #if defined(_WIN32) && !defined(__CYGWIN__) const char *homeDrive = itksys::SystemTools::GetEnv("HOMEDRIVE"); const char *homePath = itksys::SystemTools::GetEnv("HOMEPATH"); if ((homeDrive == nullptr) || (homePath == nullptr)) { itkGenericOutputMacro(<< "Environment variables HOMEDRIVE and/or HOMEPATH not set" << ". Using current working directory as home directory: " << Utf8Util::Utf8ToLocal8Bit(itksys::SystemTools::GetCurrentWorkingDirectory())); homeDirectory = Utf8Util::Utf8ToLocal8Bit(itksys::SystemTools::GetCurrentWorkingDirectory()); } else { homeDirectory = homeDrive; homeDirectory += homePath; } if (itksys::SystemTools::FileExists(Utf8Util::Local8BitToUtf8(homeDirectory).c_str()) == false) { itkGenericOutputMacro(<< "Could not find home directory at " << homeDirectory << ". Using current working directory as home directory: " << Utf8Util::Utf8ToLocal8Bit(itksys::SystemTools::GetCurrentWorkingDirectory())); homeDirectory = Utf8Util::Utf8ToLocal8Bit(itksys::SystemTools::GetCurrentWorkingDirectory()); } #else const char *home = itksys::SystemTools::GetEnv("HOME"); if (home == nullptr) { itkGenericOutputMacro(<< "Environment variable HOME not set" << ". Using current working directory as home directory: " << itksys::SystemTools::GetCurrentWorkingDirectory()); homeDirectory = itksys::SystemTools::GetCurrentWorkingDirectory(); } else homeDirectory = home; if (itksys::SystemTools::FileExists(homeDirectory.c_str()) == false) { itkGenericOutputMacro(<< "Could not find home directory at " << homeDirectory << ". Using current working directory as home directory: " << itksys::SystemTools::GetCurrentWorkingDirectory()); homeDirectory = itksys::SystemTools::GetCurrentWorkingDirectory(); } #endif // defined(_WIN32) && !defined(__CYGWIN__) optionsDirectory = homeDirectory; optionsDirectory += "/.mitk"; } optionsDirectory = itksys::SystemTools::ConvertToOutputPath(Utf8Util::Local8BitToUtf8(optionsDirectory).c_str()); if (itksys::SystemTools::CountChar(optionsDirectory.c_str(), '"') > 0) { char *unquoted = itksys::SystemTools::RemoveChars(optionsDirectory.c_str(), "\""); optionsDirectory = unquoted; delete[] unquoted; } - if (itksys::SystemTools::MakeDirectory(optionsDirectory.c_str()) == false) + if (!itksys::SystemTools::MakeDirectory(optionsDirectory.c_str()).IsSuccess()) { itkGenericExceptionMacro(<< "Could not create .mitk directory at " << Utf8Util::Utf8ToLocal8Bit(optionsDirectory)); } return Utf8Util::Utf8ToLocal8Bit(optionsDirectory); } diff --git a/Modules/Core/test/mitkBaseGeometryTest.cpp b/Modules/Core/test/mitkBaseGeometryTest.cpp index dad9c06c3a..4818e69c58 100644 --- a/Modules/Core/test/mitkBaseGeometryTest.cpp +++ b/Modules/Core/test/mitkBaseGeometryTest.cpp @@ -1,1680 +1,1680 @@ /*============================================================================ 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 "mitkTestingMacros.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class vtkMatrix4x4; class vtkMatrixToLinearTransform; class vtkLinearTransform; typedef itk::BoundingBox BoundingBox; typedef itk::BoundingBox BoundingBoxType; typedef BoundingBoxType::BoundsArrayType BoundsArrayType; typedef BoundingBoxType::Pointer BoundingBoxPointer; // Dummy instance of abstract base class class DummyTestClass : public mitk::BaseGeometry { public: DummyTestClass(){}; DummyTestClass(const DummyTestClass &other) : BaseGeometry(other){}; ~DummyTestClass() override{}; mitkClassMacro(DummyTestClass, mitk::BaseGeometry); itkNewMacro(Self); mitkNewMacro1Param(Self, const Self &); itk::LightObject::Pointer InternalClone() const override { Self::Pointer newGeometry = new Self(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } protected: void PrintSelf(std::ostream & /*os*/, itk::Indent /*indent*/) const override{}; //##Documentation //## @brief Pre- and Post-functions are empty in BaseGeometry //## //## These virtual functions allow for a different beahiour in subclasses. //## Do implement them in every subclass of BaseGeometry. If not needed, use {}. //## If this class is inherited from a subclass of BaseGeometry, call {Superclass::Pre...();};, example: // SlicedGeometry3D class void PreSetSpacing(const mitk::Vector3D &/*aSpacing*/) override{}; }; class mitkBaseGeometryTestSuite : public mitk::TestFixture { // List of Tests CPPUNIT_TEST_SUITE(mitkBaseGeometryTestSuite); // Constructor MITK_TEST(TestConstructors); MITK_TEST(TestInitialize); // Set MITK_TEST(TestSetOrigin); MITK_TEST(TestSetBounds); MITK_TEST(TestSetFloatBounds); MITK_TEST(TestSetFloatBoundsDouble); MITK_TEST(TestSetFrameOfReferenceID); MITK_TEST(TestSetIndexToWorldTransform); MITK_TEST(TestSetIndexToWorldTransformWithoutChangingSpacing); MITK_TEST(TestSetIndexToWorldTransform_WithPointerToSameTransform); MITK_TEST(TestSetSpacing); MITK_TEST(TestTransferItkToVtkTransform); MITK_TEST(TestSetIndexToWorldTransformByVtkMatrix); MITK_TEST(TestSetIdentity); MITK_TEST(TestSetImageGeometry); // Equal MITK_TEST(Equal_CloneAndOriginal_ReturnsTrue); MITK_TEST(Equal_DifferentOrigin_ReturnsFalse); MITK_TEST(Equal_DifferentIndexToWorldTransform_ReturnsFalse); MITK_TEST(Equal_DifferentSpacing_ReturnsFalse); MITK_TEST(Equal_InputIsNull_ReturnsFalse); MITK_TEST(Equal_DifferentBoundingBox_ReturnsFalse); MITK_TEST(Equal_Transforms_MinorDifferences_And_Eps); // other Functions MITK_TEST(TestComposeTransform); MITK_TEST(TestComposeVtkMatrix); MITK_TEST(TestTranslate); MITK_TEST(TestIndexToWorld); MITK_TEST(TestExecuteOperation); MITK_TEST(TestCalculateBoundingBoxRelToTransform); // MITK_TEST(TestSetTimeBounds); MITK_TEST(TestIs2DConvertable); MITK_TEST(TestGetCornerPoint); MITK_TEST(TestExtentInMM); MITK_TEST(TestGetAxisVector); MITK_TEST(TestGetCenter); MITK_TEST(TestGetDiagonalLength); MITK_TEST(TestGetExtent); MITK_TEST(TestIsInside); MITK_TEST(TestGetMatrixColumn); // test IsSubGeometry MITK_TEST(IsSubGeometry_Spacing); MITK_TEST(IsSubGeometry_TransformMatrix); MITK_TEST(IsSubGeometry_Bounds_Image); MITK_TEST(IsSubGeometry_Bounds_NoneImage); MITK_TEST(IsSubGeometry_Grid_Image); MITK_TEST(IsSubGeometry_Grid_NoneImage); MITK_TEST(IsSubGeometry_Bounds_Oblique_Image); MITK_TEST(IsSubGeometry_Bounds_Oblique_NoneImage); MITK_TEST(IsSubGeometry_Grid_Oblique_Image); MITK_TEST(IsSubGeometry_Grid_Oblique_NoneImage); CPPUNIT_TEST_SUITE_END(); // Used Variables private: mitk::Point3D aPoint; float aFloatSpacing[3]; mitk::Vector3D aSpacing; mitk::AffineTransform3D::Pointer aTransform; BoundingBoxPointer aBoundingBox; mitk::AffineTransform3D::MatrixType aMatrix; mitk::Point3D anotherPoint; mitk::Vector3D anotherSpacing; BoundingBoxPointer anotherBoundingBox; BoundingBoxPointer aThirdBoundingBox; mitk::AffineTransform3D::Pointer anotherTransform; mitk::AffineTransform3D::Pointer aThirdTransform; mitk::AffineTransform3D::MatrixType anotherMatrix; mitk::AffineTransform3D::MatrixType aThirdMatrix; DummyTestClass::Pointer aDummyGeometry; DummyTestClass::Pointer anotherDummyGeometry; DummyTestClass::Pointer aDummyGeometryOblique; public: // Set up for variables void setUp() override { mitk::FillVector3D(aFloatSpacing, 1, 1, 1); mitk::FillVector3D(aSpacing, 1, 1, 1); mitk::FillVector3D(aPoint, 0, 0, 0); // Transform aTransform = mitk::AffineTransform3D::New(); aTransform->SetIdentity(); aMatrix.SetIdentity(); anotherTransform = mitk::AffineTransform3D::New(); anotherMatrix.SetIdentity(); anotherMatrix(1, 1) = 2; anotherTransform->SetMatrix(anotherMatrix); aThirdTransform = mitk::AffineTransform3D::New(); aThirdMatrix.SetIdentity(); aThirdMatrix(1, 1) = 7; aThirdTransform->SetMatrix(aThirdMatrix); // Bounding Box float bounds[6] = { 0, 1, 0, 1, 0, 1 }; mitk::BoundingBox::BoundsArrayType b; const float* input = bounds; int j = 0; for (mitk::BoundingBox::BoundsArrayType::Iterator it = b.Begin(); j < 6; ++j) *it++ = (mitk::ScalarType) * input++; aBoundingBox = BoundingBoxType::New(); BoundingBoxType::PointsContainer::Pointer pointscontainer = BoundingBoxType::PointsContainer::New(); BoundingBoxType::PointType p; BoundingBoxType::PointIdentifier pointid; for (pointid = 0; pointid < 2; ++pointid) { unsigned int i; for (i = 0; i < 3; ++i) { p[i] = bounds[2 * i + pointid]; } pointscontainer->InsertElement(pointid, p); } aBoundingBox->SetPoints(pointscontainer); aBoundingBox->ComputeBoundingBox(); anotherBoundingBox = BoundingBoxType::New(); p[0] = 11; p[1] = 12; p[2] = 13; pointscontainer->InsertElement(1, p); anotherBoundingBox->SetPoints(pointscontainer); anotherBoundingBox->ComputeBoundingBox(); aThirdBoundingBox = BoundingBoxType::New(); p[0] = 22; p[1] = 23; p[2] = 24; pointscontainer->InsertElement(1, p); aThirdBoundingBox->SetPoints(pointscontainer); aThirdBoundingBox->ComputeBoundingBox(); mitk::FillVector3D(anotherPoint, 2, 3, 4); mitk::FillVector3D(anotherSpacing, 5, 6.5, 7); aDummyGeometry = DummyTestClass::New(); aDummyGeometry->Initialize(); anotherDummyGeometry = aDummyGeometry->Clone(); aDummyGeometryOblique = DummyTestClass::New(); aDummyGeometryOblique->Initialize(); auto newBounds = aDummyGeometryOblique->GetBounds(); newBounds[0] = 0; newBounds[1] = 5; newBounds[2] = 10; newBounds[3] = 20; newBounds[4] = 30; newBounds[5] = 40; aDummyGeometryOblique->SetBounds(newBounds); aDummyGeometryOblique->GetMatrixColumn(0); auto obliqueTransform = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::OutputVectorType rotationAxis(0.); rotationAxis[1] = 1.; obliqueTransform->Rotate3D(rotationAxis, 0.6); mitk::AffineTransform3D::OutputVectorType translation; translation[0] = 100.; translation[1] = -50.; translation[2] = -150.; obliqueTransform->SetTranslation(translation); aDummyGeometryOblique->SetIndexToWorldTransform(obliqueTransform); } void tearDown() override { aDummyGeometry = nullptr; anotherDummyGeometry = nullptr; } // Test functions void TestSetOrigin() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetOrigin(anotherPoint); CPPUNIT_ASSERT(mitk::Equal(anotherPoint, dummy->GetOrigin())); // undo changes, new and changed object need to be the same! dummy->SetOrigin(aPoint); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "TestSetOrigin"); } void TestSetImageGeometry() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetImageGeometry(true); CPPUNIT_ASSERT(dummy->GetImageGeometry()); // undo changes, new and changed object need to be the same! dummy->SetImageGeometry(false); CPPUNIT_ASSERT(dummy->GetImageGeometry() == false); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "TestSetImageGeometry"); } void TestSetFloatBounds() { float bounds[6] = {0, 11, 0, 12, 0, 13}; DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetFloatBounds(bounds); MITK_ASSERT_EQUAL(BoundingBox::ConstPointer(dummy->GetBoundingBox()), anotherBoundingBox, "BoundingBox equality"); // Wrong bounds, test needs to fail bounds[1] = 7; dummy->SetFloatBounds(bounds); MITK_ASSERT_NOT_EQUAL( BoundingBox::ConstPointer(dummy->GetBoundingBox()), anotherBoundingBox, "BoundingBox not equal"); // undo changes, new and changed object need to be the same! float originalBounds[6] = {0, 1, 0, 1, 0, 1}; dummy->SetFloatBounds(originalBounds); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Undo and equal"); } void TestSetBounds() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetBounds(anotherBoundingBox->GetBounds()); MITK_ASSERT_EQUAL(BoundingBox::ConstPointer(dummy->GetBoundingBox()), anotherBoundingBox, "Setting bounds"); // Test needs to fail now dummy->SetBounds(aThirdBoundingBox->GetBounds()); MITK_ASSERT_NOT_EQUAL( BoundingBox::ConstPointer(dummy->GetBoundingBox()), anotherBoundingBox, "Setting unequal bounds"); // undo changes, new and changed object need to be the same! dummy->SetBounds(aBoundingBox->GetBounds()); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Undo set bounds"); } void TestSetFloatBoundsDouble() { double bounds[6] = {0, 11, 0, 12, 0, 13}; DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetFloatBounds(bounds); MITK_ASSERT_EQUAL(BoundingBox::ConstPointer(dummy->GetBoundingBox()), anotherBoundingBox, "Float bounds"); // Test needs to fail now bounds[3] = 7; dummy->SetFloatBounds(bounds); MITK_ASSERT_NOT_EQUAL( BoundingBox::ConstPointer(dummy->GetBoundingBox()), anotherBoundingBox, "Float bounds unequal"); // undo changes, new and changed object need to be the same! double originalBounds[6] = {0, 1, 0, 1, 0, 1}; dummy->SetFloatBounds(originalBounds); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Undo set float bounds"); } void TestSetFrameOfReferenceID() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetFrameOfReferenceID(5); CPPUNIT_ASSERT(dummy->GetFrameOfReferenceID() == 5); // undo changes, new and changed object need to be the same! dummy->SetFrameOfReferenceID(0); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Undo set frame of reference"); } void TestSetIndexToWorldTransform() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); MITK_ASSERT_EQUAL(anotherTransform, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Compare IndexToWorldTransform 1"); // Test needs to fail now dummy->SetIndexToWorldTransform(aThirdTransform); MITK_ASSERT_NOT_EQUAL(anotherTransform, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Compare IndexToWorldTransform 2"); // undo changes, new and changed object need to be the same! dummy->SetIndexToWorldTransform(aTransform); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Compare IndexToWorldTransform 3"); } void TestSetIndexToWorldTransformWithoutChangingSpacing() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransformWithoutChangingSpacing(anotherTransform); CPPUNIT_ASSERT(mitk::Equal(aSpacing, dummy->GetSpacing(), mitk::eps, true)); // calculate a new version of anotherTransform, so that the spacing should be the same as the original spacing of // aTransform. mitk::AffineTransform3D::MatrixType::InternalMatrixType vnlmatrix; vnlmatrix = anotherTransform->GetMatrix().GetVnlMatrix(); mitk::VnlVector col; col = vnlmatrix.get_column(0).as_ref(); col.normalize(); col *= aSpacing[0]; vnlmatrix.set_column(0, col); col = vnlmatrix.get_column(1).as_ref(); col.normalize(); col *= aSpacing[1]; vnlmatrix.set_column(1, col); col = vnlmatrix.get_column(2).as_ref(); col.normalize(); col *= aSpacing[2]; vnlmatrix.set_column(2, col); mitk::Matrix3D matrix; matrix = vnlmatrix; anotherTransform->SetMatrix(matrix); CPPUNIT_ASSERT(mitk::Equal(*anotherTransform, *(dummy->GetIndexToWorldTransform()), mitk::eps, true)); } void TestSetIndexToWorldTransform_WithPointerToSameTransform() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetOrigin(anotherPoint); dummy->SetIndexToWorldTransform(anotherTransform); dummy->SetSpacing(anotherSpacing); mitk::AffineTransform3D::Pointer testTransfrom = dummy->GetIndexToWorldTransform(); mitk::Vector3D modifiedPoint = anotherPoint.GetVectorFromOrigin() * 2.; testTransfrom->SetOffset(modifiedPoint); dummy->SetIndexToWorldTransform(testTransfrom); CPPUNIT_ASSERT(mitk::Equal(modifiedPoint, dummy->GetOrigin().GetVectorFromOrigin())); } void TestSetIndexToWorldTransformByVtkMatrix() { vtkMatrix4x4 *vtkmatrix; vtkmatrix = vtkMatrix4x4::New(); vtkmatrix->Identity(); vtkmatrix->SetElement(1, 1, 2); DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransformByVtkMatrix(vtkmatrix); MITK_ASSERT_EQUAL(anotherTransform, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Compare IndexToWorldTransformByVtkMatrix 1"); // test needs to fail now vtkmatrix->SetElement(1, 1, 7); dummy->SetIndexToWorldTransformByVtkMatrix(vtkmatrix); MITK_ASSERT_NOT_EQUAL(anotherTransform, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Compare IndexToWorldTransformByVtkMatrix 2"); // undo changes, new and changed object need to be the same! vtkmatrix->SetElement(1, 1, 1); dummy->SetIndexToWorldTransformByVtkMatrix(vtkmatrix); vtkmatrix->Delete(); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Compare IndexToWorldTransformByVtkMatrix 3"); } void TestSetIdentity() { DummyTestClass::Pointer dummy = DummyTestClass::New(); // Change IndextoWorldTransform and Origin dummy->SetIndexToWorldTransform(anotherTransform); dummy->SetOrigin(anotherPoint); // Set Identity should reset ITWT and Origin dummy->SetIdentity(); MITK_ASSERT_EQUAL( aTransform, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Test set identity 1"); CPPUNIT_ASSERT(mitk::Equal(aPoint, dummy->GetOrigin())); CPPUNIT_ASSERT(mitk::Equal(aSpacing, dummy->GetSpacing())); // new and changed object need to be the same! DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Test set identity 2"); } void TestSetSpacing() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetSpacing(anotherSpacing); CPPUNIT_ASSERT(mitk::Equal(anotherSpacing, dummy->GetSpacing())); // undo changes, new and changed object need to be the same! dummy->SetSpacing(aSpacing); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Dummy spacing"); } void TestTransferItkToVtkTransform() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); // calls TransferItkToVtkTransform mitk::AffineTransform3D::Pointer dummyTransform = dummy->GetIndexToWorldTransform(); CPPUNIT_ASSERT(mitk::MatrixEqualElementWise(anotherMatrix, dummyTransform->GetMatrix())); } void TestConstructors() { // test standard constructor DummyTestClass::Pointer dummy1 = DummyTestClass::New(); bool test = dummy1->IsValid(); CPPUNIT_ASSERT(test == true); CPPUNIT_ASSERT(dummy1->GetFrameOfReferenceID() == 0); CPPUNIT_ASSERT(dummy1->GetIndexToWorldTransformLastModified() == 0); CPPUNIT_ASSERT(mitk::Equal(dummy1->GetSpacing(), aSpacing)); CPPUNIT_ASSERT(mitk::Equal(dummy1->GetOrigin(), aPoint)); CPPUNIT_ASSERT(dummy1->GetImageGeometry() == false); MITK_ASSERT_EQUAL( mitk::AffineTransform3D::Pointer(dummy1->GetIndexToWorldTransform()), aTransform, "Constructor test 1"); MITK_ASSERT_EQUAL( mitk::BaseGeometry::BoundingBoxType::ConstPointer(dummy1->GetBoundingBox()), aBoundingBox, "Constructor test 2"); DummyTestClass::Pointer dummy2 = DummyTestClass::New(); dummy2->SetOrigin(anotherPoint); float bounds[6] = {0, 11, 0, 12, 0, 13}; dummy2->SetFloatBounds(bounds); dummy2->SetIndexToWorldTransform(anotherTransform); dummy2->SetSpacing(anotherSpacing); DummyTestClass::Pointer dummy3 = DummyTestClass::New(*dummy2); MITK_ASSERT_EQUAL(dummy3, dummy2, "Dummy constructor"); } // Equal Tests void Equal_CloneAndOriginal_ReturnsTrue() { MITK_ASSERT_EQUAL(aDummyGeometry, anotherDummyGeometry, "Clone test"); } void Equal_DifferentOrigin_ReturnsFalse() { anotherDummyGeometry->SetOrigin(anotherPoint); MITK_ASSERT_NOT_EQUAL(aDummyGeometry, anotherDummyGeometry, "Different origin test"); } void Equal_DifferentIndexToWorldTransform_ReturnsFalse() { anotherDummyGeometry->SetIndexToWorldTransform(anotherTransform); MITK_ASSERT_NOT_EQUAL(aDummyGeometry, anotherDummyGeometry, "Different index to world"); } void Equal_DifferentSpacing_ReturnsFalse() { anotherDummyGeometry->SetSpacing(anotherSpacing); MITK_ASSERT_NOT_EQUAL(aDummyGeometry, anotherDummyGeometry, "Different spacing"); } void Equal_InputIsNull_ReturnsFalse() { DummyTestClass::Pointer geometryNull = nullptr; CPPUNIT_ASSERT_THROW(MITK_ASSERT_EQUAL(geometryNull, anotherDummyGeometry, "Input is null"), mitk::Exception); } void Equal_DifferentBoundingBox_ReturnsFalse() { // create different bounds to make the comparison false mitk::ScalarType bounds[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; anotherDummyGeometry->SetBounds(bounds); MITK_ASSERT_NOT_EQUAL(aDummyGeometry, anotherDummyGeometry, "Different bounding box"); } void Equal_Transforms_MinorDifferences_And_Eps() { // Verifies that the eps parameter is evaluated properly // when comparing two mitk::BaseGeometry::TransformTypes aMatrix.SetIdentity(); anotherMatrix.SetIdentity(); aMatrix(0, 1) = 0.0002; aTransform->SetMatrix(aMatrix); anotherMatrix(0, 1) = 0.0002; anotherTransform->SetMatrix(anotherMatrix); anotherTransform->SetMatrix(aMatrix); CPPUNIT_ASSERT_MESSAGE("Exact same transforms are mitk::Equal() for eps=mitk::eps", mitk::Equal(*aTransform, *anotherTransform, mitk::eps, true)); CPPUNIT_ASSERT_MESSAGE("Exact same transforms are mitk::Equal() for eps=vnl_math::eps", mitk::Equal(*aTransform, *anotherTransform, vnl_math::eps, true)); anotherMatrix(0, 1) = 0.0002 + mitk::eps; anotherTransform->SetMatrix(anotherMatrix); CPPUNIT_ASSERT_MESSAGE("Transforms of diff mitk::eps are !mitk::Equal() for eps=vnl_math::eps", !mitk::Equal(*aTransform, *anotherTransform, vnl_math::eps, true)); CPPUNIT_ASSERT_MESSAGE("Transforms of diff mitk::eps are !mitk::Equal() for eps=mitk::eps-1%", !mitk::Equal(*aTransform, *anotherTransform, mitk::eps * 0.99, true)); CPPUNIT_ASSERT_MESSAGE("Transforms of diff mitk::eps _are_ mitk::Equal() for eps=mitk::eps+1%", mitk::Equal(*aTransform, *anotherTransform, mitk::eps * 1.01, true)); } void TestComposeTransform() { // Create Transformations to set and compare mitk::AffineTransform3D::Pointer transform1; transform1 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix1; matrix1.SetIdentity(); matrix1(1, 1) = 2; transform1->SetMatrix(matrix1); // Spacing = 2 mitk::AffineTransform3D::Pointer transform2; transform2 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix2; matrix2.SetIdentity(); matrix2(1, 1) = 2; transform2->SetMatrix(matrix2); // Spacing = 2 mitk::AffineTransform3D::Pointer transform3; transform3 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix3; matrix3.SetIdentity(); matrix3(1, 1) = 4; transform3->SetMatrix(matrix3); // Spacing = 4 mitk::AffineTransform3D::Pointer transform4; transform4 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix4; matrix4.SetIdentity(); matrix4(1, 1) = 0.25; transform4->SetMatrix(matrix4); // Spacing = 0.25 // Vector to compare spacing mitk::Vector3D expectedSpacing; expectedSpacing.Fill(1.0); expectedSpacing[1] = 4; DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(transform1); // Spacing = 2 dummy->Compose(transform2); // Spacing = 4 CPPUNIT_ASSERT(mitk::Equal(dummy->GetSpacing(), expectedSpacing)); MITK_ASSERT_EQUAL( transform3, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Compose transform 2"); // 4=4 // undo changes, new and changed object need to be the same! dummy->Compose(transform4); // Spacing = 1 DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Compose transform 3"); // 1=1 } void TestComposeVtkMatrix() { // Create Transformations to set and compare mitk::AffineTransform3D::Pointer transform1; transform1 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix1; matrix1.SetIdentity(); matrix1(1, 1) = 2; transform1->SetMatrix(matrix1); // Spacing = 2 vtkMatrix4x4 *vtkmatrix2; vtkmatrix2 = vtkMatrix4x4::New(); vtkmatrix2->Identity(); vtkmatrix2->SetElement(1, 1, 2); // Spacing = 2 mitk::AffineTransform3D::Pointer transform3; transform3 = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType matrix3; matrix3.SetIdentity(); matrix3(1, 1) = 4; transform3->SetMatrix(matrix3); // Spacing = 4 vtkMatrix4x4 *vtkmatrix4; vtkmatrix4 = vtkMatrix4x4::New(); vtkmatrix4->Identity(); vtkmatrix4->SetElement(1, 1, 0.25); // Spacing = 0.25 // Vector to compare spacing mitk::Vector3D expectedSpacing; expectedSpacing.Fill(1.0); expectedSpacing[1] = 4; DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(transform1); // Spacing = 2 dummy->Compose(vtkmatrix2); // Spacing = 4 vtkmatrix2->Delete(); MITK_ASSERT_EQUAL( transform3, mitk::AffineTransform3D::Pointer(dummy->GetIndexToWorldTransform()), "Compose vtk matrix"); // 4=4 CPPUNIT_ASSERT(mitk::Equal(dummy->GetSpacing(), expectedSpacing)); // undo changes, new and changed object need to be the same! dummy->Compose(vtkmatrix4); // Spacing = 1 vtkmatrix4->Delete(); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Compose vtk"); // 1=1 } void TestTranslate() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetOrigin(anotherPoint); CPPUNIT_ASSERT(mitk::Equal(anotherPoint, dummy->GetOrigin())); // use some random values for translation mitk::Vector3D translationVector; translationVector.SetElement(0, 17.5f); translationVector.SetElement(1, -32.3f); translationVector.SetElement(2, 4.0f); // compute ground truth mitk::Point3D tmpResult = anotherPoint + translationVector; dummy->Translate(translationVector); CPPUNIT_ASSERT(mitk::Equal(dummy->GetOrigin(), tmpResult)); // undo changes translationVector *= -1; dummy->Translate(translationVector); CPPUNIT_ASSERT(mitk::Equal(dummy->GetOrigin(), anotherPoint)); // undo changes, new and changed object need to be the same! translationVector.SetElement(0, -1 * anotherPoint[0]); translationVector.SetElement(1, -1 * anotherPoint[1]); translationVector.SetElement(2, -1 * anotherPoint[2]); dummy->Translate(translationVector); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Translate test"); } // a part of the test requires axis-parallel coordinates int testIndexAndWorldConsistency(DummyTestClass::Pointer dummyGeometry) { // Testing consistency of index and world coordinate systems mitk::Point3D origin = dummyGeometry->GetOrigin(); mitk::Point3D dummyPoint; // Testing index->world->index conversion consistency dummyGeometry->WorldToIndex(origin, dummyPoint); dummyGeometry->IndexToWorld(dummyPoint, dummyPoint); CPPUNIT_ASSERT(mitk::EqualArray(dummyPoint, origin, 3, mitk::eps, true)); // Testing WorldToIndex(origin, mitk::Point3D)==(0,0,0) mitk::Point3D globalOrigin; mitk::FillVector3D(globalOrigin, 0, 0, 0); mitk::Point3D originContinuousIndex; dummyGeometry->WorldToIndex(origin, originContinuousIndex); CPPUNIT_ASSERT(mitk::EqualArray(originContinuousIndex, globalOrigin, 3, mitk::eps, true)); // Testing WorldToIndex(origin, itk::Index)==(0,0,0) itk::Index<3> itkindex; dummyGeometry->WorldToIndex(origin, itkindex); itk::Index<3> globalOriginIndex; mitk::vtk2itk(globalOrigin, globalOriginIndex); CPPUNIT_ASSERT(mitk::EqualArray(itkindex, globalOriginIndex, 3, mitk::eps, true)); // Testing WorldToIndex(origin-0.5*spacing, itk::Index)==(0,0,0) mitk::Vector3D halfSpacingStep = dummyGeometry->GetSpacing() * 0.5; mitk::Matrix3D rotation; mitk::Point3D originOffCenter = origin - halfSpacingStep; dummyGeometry->WorldToIndex(originOffCenter, itkindex); CPPUNIT_ASSERT(mitk::EqualArray(itkindex, globalOriginIndex, 3, mitk::eps, true)); // Testing WorldToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0) originOffCenter = origin + halfSpacingStep; - originOffCenter -= 0.0001; + originOffCenter -= mitk::Vector(0.0001); dummyGeometry->WorldToIndex(originOffCenter, itkindex); CPPUNIT_ASSERT(mitk::EqualArray(itkindex, globalOriginIndex, 3, mitk::eps, true)); // Testing WorldToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)"); originOffCenter = origin + halfSpacingStep; itk::Index<3> global111; mitk::FillVector3D(global111, 1, 1, 1); dummyGeometry->WorldToIndex(originOffCenter, itkindex); CPPUNIT_ASSERT(mitk::EqualArray(itkindex, global111, 3, mitk::eps, true)); // Testing WorldToIndex(GetCenter())==BoundingBox.GetCenter mitk::Point3D center = dummyGeometry->GetCenter(); mitk::Point3D centerContIndex; dummyGeometry->WorldToIndex(center, centerContIndex); mitk::BoundingBox::ConstPointer boundingBox = dummyGeometry->GetBoundingBox(); mitk::BoundingBox::PointType centerBounds = boundingBox->GetCenter(); CPPUNIT_ASSERT(mitk::Equal(centerContIndex, centerBounds)); // Testing GetCenter()==IndexToWorld(BoundingBox.GetCenter) center = dummyGeometry->GetCenter(); mitk::Point3D centerBoundsInWorldCoords; dummyGeometry->IndexToWorld(centerBounds, centerBoundsInWorldCoords); CPPUNIT_ASSERT(mitk::Equal(center, centerBoundsInWorldCoords)); // Test using random point, // Testing consistency of index and world coordinate systems mitk::Point3D point; mitk::FillVector3D(point, 3.5, -2, 4.6); // Testing index->world->index conversion consistency dummyGeometry->WorldToIndex(point, dummyPoint); dummyGeometry->IndexToWorld(dummyPoint, dummyPoint); CPPUNIT_ASSERT(mitk::EqualArray(dummyPoint, point, 3, mitk::eps, true)); return EXIT_SUCCESS; } int testIndexAndWorldConsistencyForVectors(DummyTestClass::Pointer dummyGeometry) { // Testing consistency of index and world coordinate systems for vectors mitk::Vector3D xAxisMM = dummyGeometry->GetAxisVector(0); mitk::Vector3D xAxisContinuousIndex; mitk::Point3D p, pIndex, origin; origin = dummyGeometry->GetOrigin(); p[0] = xAxisMM[0] + origin[0]; p[1] = xAxisMM[1] + origin[1]; p[2] = xAxisMM[2] + origin[2]; dummyGeometry->WorldToIndex(p, pIndex); dummyGeometry->WorldToIndex(xAxisMM, xAxisContinuousIndex); CPPUNIT_ASSERT(mitk::Equal(xAxisContinuousIndex[0], pIndex[0])); CPPUNIT_ASSERT(mitk::Equal(xAxisContinuousIndex[1], pIndex[1])); CPPUNIT_ASSERT(mitk::Equal(xAxisContinuousIndex[2], pIndex[2])); dummyGeometry->IndexToWorld(xAxisContinuousIndex, xAxisContinuousIndex); dummyGeometry->IndexToWorld(pIndex, p); CPPUNIT_ASSERT(xAxisContinuousIndex == xAxisMM); CPPUNIT_ASSERT(mitk::Equal(xAxisContinuousIndex[0], p[0] - origin[0])); CPPUNIT_ASSERT(mitk::Equal(xAxisContinuousIndex[1], p[1] - origin[1])); CPPUNIT_ASSERT(mitk::Equal(xAxisContinuousIndex[2], p[2] - origin[2])); // Test consictency for random vector mitk::Vector3D vector; mitk::FillVector3D(vector, 2.5, -3.2, 8.1); mitk::Vector3D vectorContinuousIndex; p[0] = vector[0] + origin[0]; p[1] = vector[1] + origin[1]; p[2] = vector[2] + origin[2]; dummyGeometry->WorldToIndex(p, pIndex); dummyGeometry->WorldToIndex(vector, vectorContinuousIndex); CPPUNIT_ASSERT(mitk::Equal(vectorContinuousIndex[0], pIndex[0])); CPPUNIT_ASSERT(mitk::Equal(vectorContinuousIndex[1], pIndex[1])); CPPUNIT_ASSERT(mitk::Equal(vectorContinuousIndex[2], pIndex[2])); dummyGeometry->IndexToWorld(vectorContinuousIndex, vectorContinuousIndex); dummyGeometry->IndexToWorld(pIndex, p); CPPUNIT_ASSERT(vectorContinuousIndex == vector); CPPUNIT_ASSERT(mitk::Equal(vectorContinuousIndex[0], p[0] - origin[0])); CPPUNIT_ASSERT(mitk::Equal(vectorContinuousIndex[1], p[1] - origin[1])); CPPUNIT_ASSERT(mitk::Equal(vectorContinuousIndex[2], p[2] - origin[2])); return EXIT_SUCCESS; } int testIndexAndWorldConsistencyForIndex(DummyTestClass::Pointer dummyGeometry) { // Testing consistency of index and world coordinate systems // creating testing data itk::Index<4> itkIndex4, itkIndex4b; itk::Index<3> itkIndex3, itkIndex3b; itk::Index<2> itkIndex2, itkIndex2b; itk::Index<3> mitkIndex, mitkIndexb; itkIndex4[0] = itkIndex4[1] = itkIndex4[2] = itkIndex4[3] = 4; itkIndex3[0] = itkIndex3[1] = itkIndex3[2] = 6; itkIndex2[0] = itkIndex2[1] = 2; mitkIndex[0] = mitkIndex[1] = mitkIndex[2] = 13; // check for consistency mitk::Point3D point; dummyGeometry->IndexToWorld(itkIndex2, point); dummyGeometry->WorldToIndex(point, itkIndex2b); CPPUNIT_ASSERT(((itkIndex2b[0] == itkIndex2[0]) && (itkIndex2b[1] == itkIndex2[1]))); // Testing itk::index<2> for IndexToWorld/WorldToIndex consistency dummyGeometry->IndexToWorld(itkIndex3, point); dummyGeometry->WorldToIndex(point, itkIndex3b); CPPUNIT_ASSERT( ((itkIndex3b[0] == itkIndex3[0]) && (itkIndex3b[1] == itkIndex3[1]) && (itkIndex3b[2] == itkIndex3[2]))); // Testing itk::index<3> for IndexToWorld/WorldToIndex consistency dummyGeometry->IndexToWorld(itkIndex4, point); dummyGeometry->WorldToIndex(point, itkIndex4b); CPPUNIT_ASSERT(((itkIndex4b[0] == itkIndex4[0]) && (itkIndex4b[1] == itkIndex4[1]) && (itkIndex4b[2] == itkIndex4[2]) && (itkIndex4b[3] == 0))); // Testing itk::index<3> for IndexToWorld/WorldToIndex consistency dummyGeometry->IndexToWorld(mitkIndex, point); dummyGeometry->WorldToIndex(point, mitkIndexb); CPPUNIT_ASSERT( ((mitkIndexb[0] == mitkIndex[0]) && (mitkIndexb[1] == mitkIndex[1]) && (mitkIndexb[2] == mitkIndex[2]))); // Testing mitk::Index for IndexToWorld/WorldToIndex consistency return EXIT_SUCCESS; } void TestIndexToWorld() { DummyTestClass::Pointer dummy = DummyTestClass::New(); testIndexAndWorldConsistency(dummy); testIndexAndWorldConsistencyForVectors(dummy); testIndexAndWorldConsistencyForIndex(dummy); // Geometry must not have changed DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Dummy index to world"); // Test with other geometries dummy->SetOrigin(anotherPoint); testIndexAndWorldConsistency(dummy); testIndexAndWorldConsistencyForVectors(dummy); testIndexAndWorldConsistencyForIndex(dummy); dummy->SetIndexToWorldTransform(anotherTransform); testIndexAndWorldConsistency(dummy); testIndexAndWorldConsistencyForVectors(dummy); testIndexAndWorldConsistencyForIndex(dummy); dummy->SetOrigin(anotherPoint); testIndexAndWorldConsistency(dummy); testIndexAndWorldConsistencyForVectors(dummy); testIndexAndWorldConsistencyForIndex(dummy); dummy->SetSpacing(anotherSpacing); testIndexAndWorldConsistency(dummy); testIndexAndWorldConsistencyForVectors(dummy); testIndexAndWorldConsistencyForIndex(dummy); } void TestExecuteOperation() { DummyTestClass::Pointer dummy = DummyTestClass::New(); // Do same Operations with new Dummy and compare DummyTestClass::Pointer newDummy = DummyTestClass::New(); // Test operation Nothing auto opN = new mitk::Operation(mitk::OpNOTHING); dummy->ExecuteOperation(opN); MITK_ASSERT_EQUAL(dummy, newDummy, "Dummy execute operation 1"); // Test operation Move auto opP = new mitk::PointOperation(mitk::OpMOVE, anotherPoint); dummy->ExecuteOperation(opP); CPPUNIT_ASSERT(mitk::Equal(anotherPoint, dummy->GetOrigin())); newDummy->SetOrigin(anotherPoint); MITK_ASSERT_EQUAL(dummy, newDummy, "Dummy execute operation 2"); // Test operation Scale, Scale sets spacing to scale+1 mitk::Point3D spacing; spacing[0] = anotherSpacing[0] - 1.; spacing[1] = anotherSpacing[1] - 1.; spacing[2] = anotherSpacing[2] - 1.; auto opS = new mitk::ScaleOperation(mitk::OpSCALE, spacing, anotherPoint); dummy->ExecuteOperation(opS); CPPUNIT_ASSERT(mitk::Equal(anotherSpacing, dummy->GetSpacing())); newDummy->SetSpacing(anotherSpacing); MITK_ASSERT_EQUAL(dummy, newDummy, "Dummy execute operation 3"); // change Geometry to test more cases dummy->SetIndexToWorldTransform(anotherTransform); dummy->SetSpacing(anotherSpacing); // Testing a rotation of the geometry double angle = 35.0; mitk::Vector3D rotationVector; mitk::FillVector3D(rotationVector, 1, 0, 0); mitk::Point3D center = dummy->GetCenter(); auto opR = new mitk::RotationOperation(mitk::OpROTATE, center, rotationVector, angle); dummy->ExecuteOperation(opR); mitk::Matrix3D rotation; mitk::GetRotation(dummy, rotation); mitk::Vector3D voxelStep = rotation * anotherSpacing; mitk::Vector3D voxelStepIndex; dummy->WorldToIndex(voxelStep, voxelStepIndex); mitk::Vector3D expectedVoxelStepIndex; expectedVoxelStepIndex.Fill(1); CPPUNIT_ASSERT(mitk::Equal(voxelStepIndex, expectedVoxelStepIndex)); delete opR; delete opN; delete opS; delete opP; } void TestCalculateBoundingBoxRelToTransform() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetExtentInMM(0, 15); dummy->SetExtentInMM(1, 20); dummy->SetExtentInMM(2, 8); mitk::BoundingBox::Pointer dummyBoundingBox = dummy->CalculateBoundingBoxRelativeToTransform(anotherTransform); mitk::BoundingBox::PointsContainer::Pointer pointscontainer = mitk::BoundingBox::PointsContainer::New(); mitk::BoundingBox::PointIdentifier pointid = 0; unsigned char i; mitk::AffineTransform3D::Pointer inverse = mitk::AffineTransform3D::New(); anotherTransform->GetInverse(inverse); for (i = 0; i < 8; ++i) pointscontainer->InsertElement(pointid++, inverse->TransformPoint(dummy->GetCornerPoint(i))); mitk::BoundingBox::Pointer result = mitk::BoundingBox::New(); result->SetPoints(pointscontainer); result->ComputeBoundingBox(); MITK_ASSERT_EQUAL(result, dummyBoundingBox, "BBox rel to transform"); // dummy still needs to be unchanged, except for extend DummyTestClass::Pointer newDummy = DummyTestClass::New(); newDummy->SetExtentInMM(0, 15); newDummy->SetExtentInMM(1, 20); newDummy->SetExtentInMM(2, 8); MITK_ASSERT_EQUAL(dummy, newDummy, "Dummy BBox"); } // void TestSetTimeBounds(){ // mitk::TimeBounds timeBounds; // timeBounds[0] = 1; // timeBounds[1] = 9; // DummyTestClass::Pointer dummy = DummyTestClass::New(); // dummy->SetTimeBounds(timeBounds); // mitk::TimeBounds timeBounds2 = dummy->GetTimeBounds(); // CPPUNIT_ASSERT(timeBounds[0]==timeBounds2[0]); // CPPUNIT_ASSERT(timeBounds[1]==timeBounds2[1]); // //undo changes, new and changed object need to be the same! // timeBounds[0]=mitk::ScalarTypeNumericTraits::NonpositiveMin(); // timeBounds[1]=mitk::ScalarTypeNumericTraits::max(); // DummyTestClass::Pointer newDummy = DummyTestClass::New(); // CPPUNIT_ASSERT(mitk::Equal(dummy,newDummy,mitk::eps,true)); //} void TestIs2DConvertable() { DummyTestClass::Pointer dummy = DummyTestClass::New(); // new initialized geometry is 2D convertable CPPUNIT_ASSERT(dummy->Is2DConvertable()); // Wrong Spacing needs to fail dummy->SetSpacing(anotherSpacing); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); // undo dummy->SetSpacing(aSpacing); CPPUNIT_ASSERT(dummy->Is2DConvertable()); // Wrong Origin needs to fail dummy->SetOrigin(anotherPoint); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); // undo dummy->SetOrigin(aPoint); CPPUNIT_ASSERT(dummy->Is2DConvertable()); // third dimension must not be transformed mitk::AffineTransform3D::Pointer dummyTransform = mitk::AffineTransform3D::New(); mitk::AffineTransform3D::MatrixType dummyMatrix; dummyMatrix.SetIdentity(); dummyTransform->SetMatrix(dummyMatrix); dummy->SetIndexToWorldTransform(dummyTransform); // identity matrix is 2DConvertable CPPUNIT_ASSERT(dummy->Is2DConvertable()); dummyMatrix(0, 2) = 3; dummyTransform->SetMatrix(dummyMatrix); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); dummyMatrix.SetIdentity(); dummyMatrix(1, 2) = 0.4; dummyTransform->SetMatrix(dummyMatrix); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); dummyMatrix.SetIdentity(); dummyMatrix(2, 2) = 3; dummyTransform->SetMatrix(dummyMatrix); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); dummyMatrix.SetIdentity(); dummyMatrix(2, 1) = 3; dummyTransform->SetMatrix(dummyMatrix); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); dummyMatrix.SetIdentity(); dummyMatrix(2, 0) = 3; dummyTransform->SetMatrix(dummyMatrix); CPPUNIT_ASSERT(dummy->Is2DConvertable() == false); // undo changes, new and changed object need to be the same! dummyMatrix.SetIdentity(); dummyTransform->SetMatrix(dummyMatrix); DummyTestClass::Pointer newDummy = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy, newDummy, "Is 2D convertable"); } void TestGetCornerPoint() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); double bounds[6] = {0, 11, 0, 12, 0, 13}; dummy->SetFloatBounds(bounds); mitk::Point3D corner, refCorner; // Corner 0 mitk::FillVector3D(refCorner, bounds[0], bounds[2], bounds[4]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(0); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(true, true, true); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 1 mitk::FillVector3D(refCorner, bounds[0], bounds[2], bounds[5]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(1); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(true, true, false); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 2 mitk::FillVector3D(refCorner, bounds[0], bounds[3], bounds[4]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(2); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(true, false, true); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 3 mitk::FillVector3D(refCorner, bounds[0], bounds[3], bounds[5]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(3); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(true, false, false); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 4 mitk::FillVector3D(refCorner, bounds[1], bounds[2], bounds[4]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(4); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(false, true, true); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 5 mitk::FillVector3D(refCorner, bounds[1], bounds[2], bounds[5]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(5); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(false, true, false); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 6 mitk::FillVector3D(refCorner, bounds[1], bounds[3], bounds[4]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(6); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(false, false, true); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Corner 7 mitk::FillVector3D(refCorner, bounds[1], bounds[3], bounds[5]); refCorner = anotherTransform->TransformPoint(refCorner); corner = dummy->GetCornerPoint(7); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); corner = dummy->GetCornerPoint(false, false, false); CPPUNIT_ASSERT(mitk::Equal(refCorner, corner)); // Wrong Corner needs to fail CPPUNIT_ASSERT_THROW(dummy->GetCornerPoint(20), itk::ExceptionObject); // dummy geometry must not have changed! DummyTestClass::Pointer newDummy = DummyTestClass::New(); newDummy->SetIndexToWorldTransform(anotherTransform); newDummy->SetFloatBounds(bounds); MITK_ASSERT_EQUAL(dummy, newDummy, "Corner point"); } void TestExtentInMM() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetExtentInMM(0, 50); CPPUNIT_ASSERT(mitk::Equal(50., dummy->GetExtentInMM(0))); // Vnl Matrix has changed. The next line only works because the spacing is 1! CPPUNIT_ASSERT( mitk::Equal(50., dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0).magnitude())); // Smaller extent than original dummy->SetExtentInMM(0, 5); CPPUNIT_ASSERT(mitk::Equal(5., dummy->GetExtentInMM(0))); CPPUNIT_ASSERT( mitk::Equal(5., dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(0).magnitude())); dummy->SetExtentInMM(1, 4); CPPUNIT_ASSERT(mitk::Equal(4., dummy->GetExtentInMM(1))); CPPUNIT_ASSERT( mitk::Equal(4., dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(1).magnitude())); dummy->SetExtentInMM(2, 2.5); CPPUNIT_ASSERT(mitk::Equal(2.5, dummy->GetExtentInMM(2))); CPPUNIT_ASSERT( mitk::Equal(2.5, dummy->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix().get_column(2).magnitude())); } void TestGetAxisVector() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); double bounds[6] = {0, 11, 0, 12, 0, 13}; dummy->SetFloatBounds(bounds); mitk::Vector3D vector; mitk::FillVector3D(vector, bounds[1], 0, 0); dummy->IndexToWorld(vector, vector); CPPUNIT_ASSERT(mitk::Equal(dummy->GetAxisVector(0), vector)); mitk::FillVector3D(vector, 0, bounds[3], 0); dummy->IndexToWorld(vector, vector); CPPUNIT_ASSERT(mitk::Equal(dummy->GetAxisVector(1), vector)); mitk::FillVector3D(vector, 0, 0, bounds[5]); dummy->IndexToWorld(vector, vector); CPPUNIT_ASSERT(mitk::Equal(dummy->GetAxisVector(2), vector)); } void TestGetCenter() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); double bounds[6] = {0, 11, 2, 12, 1, 13}; dummy->SetFloatBounds(bounds); mitk::Point3D refCenter; for (int i = 0; i < 3; i++) refCenter.SetElement(i, (bounds[2 * i] + bounds[2 * i + 1]) / 2.0); dummy->IndexToWorld(refCenter, refCenter); CPPUNIT_ASSERT(mitk::Equal(dummy->GetCenter(), refCenter)); } void TestGetDiagonalLength() { DummyTestClass::Pointer dummy = DummyTestClass::New(); double bounds[6] = {1, 3, 5, 8, 7.5, 11.5}; dummy->SetFloatBounds(bounds); // 3-1=2, 8-5=3, 11.5-7.5=4; 2^2+3^2+4^2 = 29 double expectedLength = sqrt(29.); CPPUNIT_ASSERT(mitk::Equal(expectedLength, dummy->GetDiagonalLength(), mitk::eps, true)); CPPUNIT_ASSERT(mitk::Equal(29., dummy->GetDiagonalLength2(), mitk::eps, true)); // dummy must not have changed DummyTestClass::Pointer newDummy = DummyTestClass::New(); newDummy->SetFloatBounds(bounds); MITK_ASSERT_EQUAL(dummy, newDummy, "Diagonal length"); } void TestGetExtent() { DummyTestClass::Pointer dummy = DummyTestClass::New(); double bounds[6] = {1, 3, 5, 8, 7.5, 11.5}; dummy->SetFloatBounds(bounds); CPPUNIT_ASSERT(mitk::Equal(2., dummy->GetExtent(0))); CPPUNIT_ASSERT(mitk::Equal(3., dummy->GetExtent(1))); CPPUNIT_ASSERT(mitk::Equal(4., dummy->GetExtent(2))); // dummy must not have changed DummyTestClass::Pointer newDummy = DummyTestClass::New(); newDummy->SetFloatBounds(bounds); MITK_ASSERT_EQUAL(dummy, newDummy, "Extend"); } void TestIsInside() { DummyTestClass::Pointer dummy = DummyTestClass::New(); double bounds[6] = {1, 3, 5, 8, 7.5, 11.5}; dummy->SetFloatBounds(bounds); mitk::Point3D insidePoint; mitk::Point3D outsidePoint; mitk::FillVector3D(insidePoint, 2, 6, 7.6); mitk::FillVector3D(outsidePoint, 0, 9, 8.2); CPPUNIT_ASSERT(dummy->IsIndexInside(insidePoint)); CPPUNIT_ASSERT(false == dummy->IsIndexInside(outsidePoint)); dummy->IndexToWorld(insidePoint, insidePoint); dummy->IndexToWorld(outsidePoint, outsidePoint); CPPUNIT_ASSERT(dummy->IsInside(insidePoint)); CPPUNIT_ASSERT(false == dummy->IsInside(outsidePoint)); // dummy must not have changed DummyTestClass::Pointer newDummy = DummyTestClass::New(); newDummy->SetFloatBounds(bounds); MITK_ASSERT_EQUAL(dummy, newDummy, "Is inside"); } void TestInitialize() { // test standard constructor DummyTestClass::Pointer dummy1 = DummyTestClass::New(); DummyTestClass::Pointer dummy2 = DummyTestClass::New(); dummy2->SetOrigin(anotherPoint); dummy2->SetBounds(anotherBoundingBox->GetBounds()); // mitk::TimeBounds timeBounds; // timeBounds[0] = 1; // timeBounds[1] = 9; // dummy2->SetTimeBounds(timeBounds); dummy2->SetIndexToWorldTransform(anotherTransform); dummy2->SetSpacing(anotherSpacing); dummy1->InitializeGeometry(dummy2); MITK_ASSERT_EQUAL(dummy1, dummy2, "Initialize 1"); dummy1->Initialize(); DummyTestClass::Pointer dummy3 = DummyTestClass::New(); MITK_ASSERT_EQUAL(dummy3, dummy1, "Initialize 2"); } void TestGetMatrixColumn() { DummyTestClass::Pointer dummy = DummyTestClass::New(); dummy->SetIndexToWorldTransform(anotherTransform); mitk::Vector3D testVector, refVector; testVector.SetVnlVector(dummy->GetMatrixColumn(0)); mitk::FillVector3D(refVector, 1, 0, 0); CPPUNIT_ASSERT(testVector == refVector); testVector.SetVnlVector(dummy->GetMatrixColumn(1)); mitk::FillVector3D(refVector, 0, 2, 0); CPPUNIT_ASSERT(testVector == refVector); testVector.SetVnlVector(dummy->GetMatrixColumn(2)); mitk::FillVector3D(refVector, 0, 0, 1); CPPUNIT_ASSERT(testVector == refVector); // dummy must not have changed DummyTestClass::Pointer newDummy = DummyTestClass::New(); newDummy->SetIndexToWorldTransform(anotherTransform); MITK_ASSERT_EQUAL(dummy, newDummy, "GetMatrixColumn"); } void IsSubGeometry_Spacing() { CPPUNIT_ASSERT(mitk::IsSubGeometry(*aDummyGeometry, *aDummyGeometry, mitk::eps, true)); for (unsigned int i = 0; i < 3; ++i) { mitk::Vector3D wrongSpacing = aDummyGeometry->GetSpacing(); wrongSpacing[i] += mitk::eps * 2; auto wrongGeometry = aDummyGeometry->Clone(); wrongGeometry->SetSpacing(wrongSpacing); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometry, mitk::eps, true)); } for (unsigned int i = 0; i < 3; ++i) { mitk::Vector3D wrongSpacing = aDummyGeometry->GetSpacing(); wrongSpacing[i] -= mitk::eps * 2; auto wrongGeometry = aDummyGeometry->Clone(); wrongGeometry->SetSpacing(wrongSpacing); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometry, mitk::eps, true)); } } void IsSubGeometry_TransformMatrix() { CPPUNIT_ASSERT(mitk::IsSubGeometry(*aDummyGeometry, *aDummyGeometry, mitk::eps, true)); for (unsigned int i = 0; i < 3; ++i) { for (unsigned int j = 0; j < 3; ++j) { itk::Matrix wrongMatrix = aDummyGeometry->GetIndexToWorldTransform()->GetMatrix(); wrongMatrix[i][j] += mitk::eps * 2; auto wrongGeometry = aDummyGeometry->Clone(); wrongGeometry->GetIndexToWorldTransform()->SetMatrix(wrongMatrix); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometry, mitk::eps, true)); } } } void IsSubGeometry_Bounds_NoneImage() { IsSubGeometry_Bounds_internal(false); } void IsSubGeometry_Bounds_Image() { IsSubGeometry_Bounds_internal(true); } void IsSubGeometry_Bounds_internal(bool isImage) { auto newBounds = aDummyGeometry->GetBounds(); newBounds[0] = 10; newBounds[1] = 20; newBounds[2] = 10; newBounds[3] = 20; newBounds[4] = 10; newBounds[5] = 20; aDummyGeometry->SetBounds(newBounds); aDummyGeometry->SetImageGeometry(isImage); CPPUNIT_ASSERT(mitk::IsSubGeometry(*aDummyGeometry, *aDummyGeometry, mitk::eps, true)); for (unsigned int i = 0; i < 6; ++i) { auto legalBounds = newBounds; if (i % 2 == 0) { legalBounds[i] += 1; } else { legalBounds[i] -= 1; } auto legalGeometry = aDummyGeometry->Clone(); legalGeometry->SetBounds(legalBounds); CPPUNIT_ASSERT(mitk::IsSubGeometry(*legalGeometry, *aDummyGeometry, mitk::eps, true)); } for (unsigned int i = 0; i < 6; ++i) { auto wrongBounds = aDummyGeometry->GetBounds(); if (i % 2 == 0) { wrongBounds[i] -= 1; } else { wrongBounds[i] += 1; } auto wrongGeometry = aDummyGeometry->Clone(); wrongGeometry->SetBounds(wrongBounds); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometry, mitk::eps, true)); } } void IsSubGeometry_Grid_Image() { IsSubGeometry_Grid_internal(true); } void IsSubGeometry_Grid_NoneImage() { IsSubGeometry_Grid_internal(false); } void IsSubGeometry_Grid_internal(bool isImage) { auto newBounds = aDummyGeometry->GetBounds(); newBounds[0] = 0; newBounds[1] = 20; newBounds[2] = 0; newBounds[3] = 20; newBounds[4] = 0; newBounds[5] = 20; aDummyGeometry->SetBounds(newBounds); aDummyGeometry->SetImageGeometry(isImage); auto smallerGeometry = aDummyGeometry->Clone(); newBounds[0] = 5; newBounds[1] = 10; newBounds[2] = 5; newBounds[3] = 10; newBounds[4] = 5; newBounds[5] = 10; smallerGeometry->SetBounds(newBounds); //legal negative shift for (unsigned int i = 0; i < 3; ++i) { auto legalOrigin = smallerGeometry->GetOrigin(); legalOrigin[i] -= smallerGeometry->GetSpacing()[i]; auto legalGeometry = smallerGeometry->Clone(); legalGeometry->SetOrigin(legalOrigin); CPPUNIT_ASSERT(mitk::IsSubGeometry(*legalGeometry, *aDummyGeometry, mitk::eps, true)); } //legal positive shift for (unsigned int i = 0; i < 3; ++i) { auto legalOrigin = smallerGeometry->GetOrigin(); legalOrigin[i] += smallerGeometry->GetSpacing()[i]; auto legalGeometry = smallerGeometry->Clone(); legalGeometry->SetOrigin(legalOrigin); CPPUNIT_ASSERT(mitk::IsSubGeometry(*legalGeometry, *aDummyGeometry, mitk::eps, true)); } //wrong negative shift for (unsigned int i = 0; i < 3; ++i) { auto wrongOrigin = smallerGeometry->GetOrigin(); wrongOrigin[i] -= 2 * mitk::eps; auto wrongGeometry = smallerGeometry->Clone(); wrongGeometry->SetOrigin(wrongOrigin); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometry, mitk::eps, true)); } //wrong positive shift for (unsigned int i = 0; i < 3; ++i) { auto wrongOrigin = smallerGeometry->GetOrigin(); wrongOrigin[i] += 2 * mitk::eps; auto wrongGeometry = smallerGeometry->Clone(); wrongGeometry->SetOrigin(wrongOrigin); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometry, mitk::eps, true)); } } void IsSubGeometry_Bounds_Oblique_NoneImage() { IsSubGeometry_Bounds_Oblique_internal(false); } void IsSubGeometry_Bounds_Oblique_Image() { IsSubGeometry_Bounds_Oblique_internal(true); } void IsSubGeometry_Bounds_Oblique_internal(bool isImage) { auto newBounds = aDummyGeometryOblique->GetBounds(); aDummyGeometryOblique->SetImageGeometry(isImage); //REMARK: used NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION to compensate rounding errors that //are interoduced when transforming points/indeces due to the oblique geometry. CPPUNIT_ASSERT(mitk::IsSubGeometry(*aDummyGeometryOblique, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); for (unsigned int i = 0; i < 6; ++i) { auto legalBounds = newBounds; if (i % 2 == 0) { legalBounds[i] += 1; } else { legalBounds[i] -= 1; } auto legalGeometry = aDummyGeometryOblique->Clone(); legalGeometry->SetBounds(legalBounds); CPPUNIT_ASSERT(mitk::IsSubGeometry(*legalGeometry, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); } for (unsigned int i = 0; i < 6; ++i) { auto wrongBounds = newBounds; if (i % 2 == 0) { wrongBounds[i] -= 1; } else { wrongBounds[i] += 1; } auto wrongGeometry = aDummyGeometryOblique->Clone(); wrongGeometry->SetBounds(wrongBounds); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); } } void IsSubGeometry_Grid_Oblique_NoneImage() { IsSubGeometry_Grid_Oblique_internal(false); } void IsSubGeometry_Grid_Oblique_Image() { IsSubGeometry_Grid_Oblique_internal(true); } void IsSubGeometry_Grid_Oblique_internal(bool isImage) { auto newBounds = aDummyGeometryOblique->GetBounds(); newBounds[0] = 0; newBounds[1] = 20; newBounds[2] = 0; newBounds[3] = 20; newBounds[4] = 0; newBounds[5] = 20; aDummyGeometryOblique->SetBounds(newBounds); aDummyGeometryOblique->SetImageGeometry(isImage); auto smallerGeometry = aDummyGeometryOblique->Clone(); newBounds[0] = 5; newBounds[1] = 10; newBounds[2] = 5; newBounds[3] = 10; newBounds[4] = 5; newBounds[5] = 10; smallerGeometry->SetBounds(newBounds); //REMARK: used NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION in the following checks //to compensate rounding errors that are interoduced when transforming points/indeces //due to the oblique geometry. //legal negative shift for (unsigned int i = 0; i < 3; ++i) { auto legalOrigin = smallerGeometry->GetOrigin(); mitk::Point3D index; smallerGeometry->WorldToIndex(legalOrigin, index); index[i] -= 1; smallerGeometry->IndexToWorld(index, legalOrigin); auto legalGeometry = smallerGeometry->Clone(); legalGeometry->SetOrigin(legalOrigin); CPPUNIT_ASSERT(mitk::IsSubGeometry(*legalGeometry, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); } //legal positive shift for (unsigned int i = 0; i < 3; ++i) { auto legalOrigin = smallerGeometry->GetOrigin(); mitk::Point3D index; smallerGeometry->WorldToIndex(legalOrigin, index); index[i] += 1; smallerGeometry->IndexToWorld(index, legalOrigin); auto legalGeometry = smallerGeometry->Clone(); legalGeometry->SetOrigin(legalOrigin); CPPUNIT_ASSERT(mitk::IsSubGeometry(*legalGeometry, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); } //wrong negative shift for (unsigned int i = 0; i < 3; ++i) { auto wrongOrigin = smallerGeometry->GetOrigin(); wrongOrigin[i] -= 2 * mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION; auto wrongGeometry = smallerGeometry->Clone(); wrongGeometry->SetOrigin(wrongOrigin); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); } //wrong positive shift for (unsigned int i = 0; i < 3; ++i) { auto wrongOrigin = smallerGeometry->GetOrigin(); wrongOrigin[i] += 2 * mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION; auto wrongGeometry = smallerGeometry->Clone(); wrongGeometry->SetOrigin(wrongOrigin); CPPUNIT_ASSERT(!mitk::IsSubGeometry(*wrongGeometry, *aDummyGeometryOblique, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION, mitk::NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION, true)); } } }; // end class mitkBaseGeometryTestSuite MITK_TEST_SUITE_REGISTRATION(mitkBaseGeometry) diff --git a/Modules/Core/test/mitkExtractSliceFilterTest.cpp b/Modules/Core/test/mitkExtractSliceFilterTest.cpp index 3c726ccdd9..afe25f5d1f 100644 --- a/Modules/Core/test/mitkExtractSliceFilterTest.cpp +++ b/Modules/Core/test/mitkExtractSliceFilterTest.cpp @@ -1,1069 +1,1069 @@ /*============================================================================ 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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // use this to create the test volume on the fly #define CREATE_VOLUME // use this to save the created volume //#define SAVE_VOLUME // use this to calculate the error from the sphere mathematical model to our pixel based one //#define CALC_TESTFAILURE_DEVIATION // use this to render an oblique slice through a specified image //#define SHOW_SLICE_IN_RENDER_WINDOW // use this to have infos printed in mitkLog //#define EXTRACTOR_DEBUG /*these are the deviations calculated by the function CalcTestFailureDeviation (see for details)*/ #define Testfailure_Deviation_Mean_128 0.853842 #define Testfailure_Deviation_Volume_128 0.145184 #define Testfailure_Deviation_Diameter_128 1.5625 #define Testfailure_Deviation_Mean_256 0.397693 #define Testfailure_Deviation_Volume_256 0.0141357 #define Testfailure_Deviation_Diameter_256 0.78125 #define Testfailure_Deviation_Mean_512 0.205277 #define Testfailure_Deviation_Volume_512 0.01993 #define Testfailure_Deviation_Diameter_512 0.390625 class mitkExtractSliceFilterTestClass { public: static void TestSlice(mitk::PlaneGeometry *planeGeometry, std::string testname) { TestPlane = planeGeometry; TestName = testname; mitk::ScalarType centerCoordValue = TestvolumeSize / 2.0; mitk::ScalarType center[3] = {centerCoordValue, centerCoordValue, centerCoordValue}; mitk::Point3D centerIndex(center); double radius = TestvolumeSize / 4.0; if (TestPlane->Distance(centerIndex) >= radius) return; // outside sphere // feed ExtractSliceFilter mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(); slicer->SetInput(TestVolume); slicer->SetWorldGeometry(TestPlane); slicer->Update(); MITK_TEST_CONDITION_REQUIRED(slicer->GetOutput() != nullptr, "Extractor returned a slice"); mitk::Image::Pointer reslicedImage = slicer->GetOutput(); AccessFixedDimensionByItk(reslicedImage, TestSphereRadiusByItk, 2); AccessFixedDimensionByItk(reslicedImage, TestSphereAreaByItk, 2); /* double devArea, devDiameter; if(TestvolumeSize == 128.0){ devArea = Testfailure_Deviation_Volume_128; devDiameter = Testfailure_Deviation_Diameter_128; } else if(TestvolumeSize == 256.0){devArea = Testfailure_Deviation_Volume_256; devDiameter = Testfailure_Deviation_Diameter_256;} else if (TestvolumeSize == 512.0){devArea = Testfailure_Deviation_Volume_512; devDiameter = Testfailure_Deviation_Diameter_512;} else{devArea = Testfailure_Deviation_Volume_128; devDiameter = Testfailure_Deviation_Diameter_128;} */ std::string areatestName = TestName.append(" area"); std::string diametertestName = TestName.append(" testing diameter"); // TODO think about the deviation, 1% makes no sense at all MITK_TEST_CONDITION(std::abs(100 - testResults.percentageAreaCalcToPixel) < 1, areatestName); MITK_TEST_CONDITION(std::abs(100 - testResults.percentageRadiusToPixel) < 1, diametertestName); #ifdef EXTRACTOR_DEBUG MITK_INFO << TestName << " >>> " << "planeDistanceToSphereCenter: " << testResults.planeDistanceToSphereCenter; MITK_INFO << "area in pixels: " << testResults.areaInPixel << " <-> area in mm: " << testResults.areaCalculated << " = " << testResults.percentageAreaCalcToPixel << "%"; MITK_INFO << "calculated diameter: " << testResults.diameterCalculated << " <-> diameter in mm: " << testResults.diameterInMM << " <-> diameter in pixel: " << testResults.diameterInPixel << " = " << testResults.percentageRadiusToPixel << "%"; #endif } /* * get the radius of the slice of a sphere based on pixel distance from edge to edge of the circle. */ template static void TestSphereRadiusByItk(itk::Image *inputImage) { typedef itk::Image InputImageType; // set the index to the middle of the image's edge at x and y axis typename InputImageType::IndexType currentIndexX; currentIndexX[0] = (int)(TestvolumeSize / 2.0); currentIndexX[1] = 0; typename InputImageType::IndexType currentIndexY; currentIndexY[0] = 0; currentIndexY[1] = (int)(TestvolumeSize / 2.0); // remember the last pixel value double lastValueX = inputImage->GetPixel(currentIndexX); double lastValueY = inputImage->GetPixel(currentIndexY); // storage for the index marks std::vector indicesX; std::vector indicesY; /*Get four indices on the edge of the circle*/ while (currentIndexX[1] < TestvolumeSize && currentIndexX[0] < TestvolumeSize) { // move x direction currentIndexX[1] += 1; // move y direction currentIndexY[0] += 1; if (inputImage->GetPixel(currentIndexX) > lastValueX) { // mark the current index typename InputImageType::IndexType markIndex; markIndex[0] = currentIndexX[0]; markIndex[1] = currentIndexX[1]; indicesX.push_back(markIndex); } else if (inputImage->GetPixel(currentIndexX) < lastValueX) { // mark the current index typename InputImageType::IndexType markIndex; markIndex[0] = currentIndexX[0]; markIndex[1] = currentIndexX[1] - 1; // value inside the sphere indicesX.push_back(markIndex); } if (inputImage->GetPixel(currentIndexY) > lastValueY) { // mark the current index typename InputImageType::IndexType markIndex; markIndex[0] = currentIndexY[0]; markIndex[1] = currentIndexY[1]; indicesY.push_back(markIndex); } else if (inputImage->GetPixel(currentIndexY) < lastValueY) { // mark the current index typename InputImageType::IndexType markIndex; markIndex[0] = currentIndexY[0]; markIndex[1] = currentIndexY[1] - 1; // value inside the sphere indicesY.push_back(markIndex); } // found both marks? if (indicesX.size() == 2 && indicesY.size() == 2) break; // the new 'last' values lastValueX = inputImage->GetPixel(currentIndexX); lastValueY = inputImage->GetPixel(currentIndexY); } /* *If we are here we found the four marks on the edge of the circle. *For the case our plane is rotated and shifted, we have to calculate the center of the circle, *else the center is the intersection of both straight lines between the marks. *When we have the center, the diameter of the circle will be checked to the reference value(math!). */ // each distance from the first mark of each direction to the center of the straight line between the marks double distanceToCenterX = std::abs(indicesX[0][1] - indicesX[1][1]) / 2.0; // double distanceToCenterY = std::abs(indicesY[0][0] - indicesY[1][0]) / 2.0; // the center of the straight lines typename InputImageType::IndexType centerX; // typename InputImageType::IndexType centerY; centerX[0] = indicesX[0][0]; centerX[1] = indicesX[0][1] + distanceToCenterX; // TODO think about implicit cast to int. this is not the real center of the image, which could be between two // pixels // centerY[0] = indicesY[0][0] + distanceToCenterY; // centerY[1] = inidcesY[0][1]; typename InputImageType::IndexType currentIndex(centerX); lastValueX = inputImage->GetPixel(currentIndex); long sumpixels = 0; std::vector diameterIndices; // move up while (currentIndex[1] < TestvolumeSize) { currentIndex[1] += 1; if (inputImage->GetPixel(currentIndex) != lastValueX) { typename InputImageType::IndexType markIndex; markIndex[0] = currentIndex[0]; markIndex[1] = currentIndex[1] - 1; diameterIndices.push_back(markIndex); break; } sumpixels++; lastValueX = inputImage->GetPixel(currentIndex); } currentIndex[1] -= sumpixels; // move back to center to go in the other direction lastValueX = inputImage->GetPixel(currentIndex); // move down while (currentIndex[1] >= 0) { currentIndex[1] -= 1; if (inputImage->GetPixel(currentIndex) != lastValueX) { typename InputImageType::IndexType markIndex; markIndex[0] = currentIndex[0]; markIndex[1] = currentIndex[1]; // outside sphere because we want to calculate the distance from edge to edge diameterIndices.push_back(markIndex); break; } sumpixels++; lastValueX = inputImage->GetPixel(currentIndex); } /* *Now sumpixels should be the apromximate diameter of the circle. This is checked with the calculated diameter from *the plane transformation(math). */ mitk::Point3D volumeCenter; volumeCenter[0] = volumeCenter[1] = volumeCenter[2] = TestvolumeSize / 2.0; double planeDistanceToSphereCenter = TestPlane->Distance(volumeCenter); double sphereRadius = TestvolumeSize / 4.0; // calculate the radius of the circle cut from the sphere by the plane double diameter = 2.0 * std::sqrt(std::pow(sphereRadius, 2) - std::pow(planeDistanceToSphereCenter, 2)); double percentageRadiusToPixel = 100 / diameter * sumpixels; /* *calculate the radius in mm by the both marks of the center line by using the world coordinates */ // get the points as 3D coordinates mitk::Vector3D diameterPointRight, diameterPointLeft; diameterPointRight[2] = diameterPointLeft[2] = 0.0; diameterPointLeft[0] = diameterIndices[0][0]; diameterPointLeft[1] = diameterIndices[0][1]; diameterPointRight[0] = diameterIndices[1][0]; diameterPointRight[1] = diameterIndices[1][1]; // transform to worldcoordinates TestVolume->GetGeometry()->IndexToWorld(diameterPointLeft, diameterPointLeft); TestVolume->GetGeometry()->IndexToWorld(diameterPointRight, diameterPointRight); // euklidian distance double diameterInMM = ((diameterPointLeft * -1.0) + diameterPointRight).GetNorm(); testResults.diameterInMM = diameterInMM; testResults.diameterCalculated = diameter; testResults.diameterInPixel = sumpixels; testResults.percentageRadiusToPixel = percentageRadiusToPixel; testResults.planeDistanceToSphereCenter = planeDistanceToSphereCenter; } /*brute force the area pixel by pixel*/ template static void TestSphereAreaByItk(itk::Image *inputImage) { typedef itk::Image InputImageType; typedef itk::ImageRegionConstIterator ImageIterator; ImageIterator imageIterator(inputImage, inputImage->GetLargestPossibleRegion()); imageIterator.GoToBegin(); int sumPixelsInArea = 0; while (!imageIterator.IsAtEnd()) { if (inputImage->GetPixel(imageIterator.GetIndex()) == pixelValueSet) sumPixelsInArea++; ++imageIterator; } mitk::Point3D volumeCenter; volumeCenter[0] = volumeCenter[1] = volumeCenter[2] = TestvolumeSize / 2.0; double planeDistanceToSphereCenter = TestPlane->Distance(volumeCenter); double sphereRadius = TestvolumeSize / 4.0; // calculate the radius of the circle cut from the sphere by the plane double radius = std::sqrt(std::pow(sphereRadius, 2) - std::pow(planeDistanceToSphereCenter, 2)); double areaInMM = 3.14159265358979 * std::pow(radius, 2); testResults.areaCalculated = areaInMM; testResults.areaInPixel = sumPixelsInArea; testResults.percentageAreaCalcToPixel = 100 / areaInMM * sumPixelsInArea; } /* * random a voxel. define plane through this voxel. reslice at the plane. compare the pixel vaues of the voxel * in the volume with the pixel value in the resliced image. * there are some indice shifting problems which causes the test to fail for oblique planes. seems like the chosen * worldcoordinate is not corresponding to the index in the 2D image. and so the pixel values are not the same as * expected. */ static void PixelvalueBasedTest() { /* setup itk image */ typedef itk::Image ImageType; typedef itk::ImageRegionConstIterator ImageIterator; ImageType::Pointer image = ImageType::New(); ImageType::IndexType start; start[0] = start[1] = start[2] = 0; ImageType::SizeType size; size[0] = size[1] = size[2] = 32; ImageType::RegionType imgRegion; imgRegion.SetSize(size); imgRegion.SetIndex(start); image->SetRegions(imgRegion); - image->SetSpacing(1.0); + image->SetSpacing(mitk::Vector(1.0)); image->Allocate(); ImageIterator imageIterator(image, image->GetLargestPossibleRegion()); imageIterator.GoToBegin(); unsigned short pixelValue = 0; // fill the image with distinct values while (!imageIterator.IsAtEnd()) { image->SetPixel(imageIterator.GetIndex(), pixelValue); ++imageIterator; ++pixelValue; } /* end setup itk image */ mitk::Image::Pointer imageInMitk; CastToMitkImage(image, imageInMitk); /*mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New(); writer->SetInput(imageInMitk); std::string file = "C:\\Users\\schroedt\\Desktop\\cube.nrrd"; writer->SetFileName(file); writer->Update();*/ PixelvalueBasedTestByPlane(imageInMitk, mitk::AnatomicalPlane::Coronal); PixelvalueBasedTestByPlane(imageInMitk, mitk::AnatomicalPlane::Sagittal); PixelvalueBasedTestByPlane(imageInMitk, mitk::AnatomicalPlane::Axial); } static void PixelvalueBasedTestByPlane(mitk::Image *imageInMitk, mitk::AnatomicalPlane orientation) { typedef itk::Image ImageType; // set the seed of the rand function srand((unsigned)time(nullptr)); /* setup a random orthogonal plane */ int sliceindex = 17; // rand() % 32; bool isFrontside = true; bool isRotated = false; if (orientation == mitk::AnatomicalPlane::Axial) { /*isFrontside = false; isRotated = true;*/ } mitk::PlaneGeometry::Pointer plane = mitk::PlaneGeometry::New(); plane->InitializeStandardPlane(imageInMitk->GetGeometry(), orientation, sliceindex, isFrontside, isRotated); mitk::Point3D origin = plane->GetOrigin(); mitk::Vector3D normal; normal = plane->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 plane->SetOrigin(origin); // we dont need this any more, because we are only testing orthogonal planes /*mitk::Vector3D rotationVector; rotationVector[0] = randFloat(); rotationVector[1] = randFloat(); rotationVector[2] = randFloat(); float degree = randFloat() * 180.0; mitk::RotationOperation* op = new mitk::RotationOperation(mitk::OpROTATE, plane->GetCenter(), rotationVector, degree); plane->ExecuteOperation(op); delete op;*/ /* end setup plane */ /* define a point in the 3D volume. * add the two axis vectors of the plane (each multiplied with a * random number) to the origin. now the two random numbers * become our index coordinates in the 2D image, because the * length of the axis vectors is 1. */ mitk::Point3D planeOrigin = plane->GetOrigin(); mitk::Vector3D axis0, axis1; axis0 = plane->GetAxisVector(0); axis1 = plane->GetAxisVector(1); axis0.Normalize(); axis1.Normalize(); unsigned char n1 = 7; // rand() % 32; unsigned char n2 = 13; // rand() % 32; mitk::Point3D testPoint3DInWorld; testPoint3DInWorld = planeOrigin + (axis0 * n1) + (axis1 * n2); // get the index of the point in the 3D volume ImageType::IndexType testPoint3DInIndex; imageInMitk->GetGeometry()->WorldToIndex(testPoint3DInWorld, testPoint3DInIndex); itk::Index<3> testPoint2DInIndex; /* end define a point in the 3D volume.*/ // do reslicing at the plane mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(); slicer->SetInput(imageInMitk); slicer->SetWorldGeometry(plane); slicer->Update(); mitk::Image::Pointer slice = slicer->GetOutput(); // Get TestPoiont3D as Index in Slice slice->GetGeometry()->WorldToIndex(testPoint3DInWorld, testPoint2DInIndex); mitk::Point3D p, sliceIndexToWorld, imageIndexToWorld; p[0] = testPoint2DInIndex[0]; p[1] = testPoint2DInIndex[1]; p[2] = testPoint2DInIndex[2]; slice->GetGeometry()->IndexToWorld(p, sliceIndexToWorld); p[0] = testPoint3DInIndex[0]; p[1] = testPoint3DInIndex[1]; p[2] = testPoint3DInIndex[2]; imageInMitk->GetGeometry()->IndexToWorld(p, imageIndexToWorld); itk::Index<2> testPoint2DIn2DIndex; testPoint2DIn2DIndex[0] = testPoint2DInIndex[0]; testPoint2DIn2DIndex[1] = testPoint2DInIndex[1]; typedef mitk::ImagePixelReadAccessor VolumeReadAccessorType; typedef mitk::ImagePixelReadAccessor SliceReadAccessorType; VolumeReadAccessorType VolumeReadAccessor(imageInMitk); SliceReadAccessorType SliceReadAccessor(slice); // compare the pixelvalues of the defined point in the 3D volume with the value of the resliced image unsigned short valueAt3DVolume = VolumeReadAccessor.GetPixelByIndex(testPoint3DInIndex); unsigned short valueAtSlice = SliceReadAccessor.GetPixelByIndex(testPoint2DIn2DIndex); // valueAt3DVolume == valueAtSlice is not always working. because of rounding errors // indices are shifted MITK_TEST_CONDITION(valueAt3DVolume == valueAtSlice, "comparing pixelvalues for orthogonal plane"); vtkSmartPointer imageInVtk = imageInMitk->GetVtkImageData(); vtkSmartPointer sliceInVtk = slice->GetVtkImageData(); double PixelvalueByMitkOutput = sliceInVtk->GetScalarComponentAsDouble(n1, n2, 0, 0); // double valueVTKinImage = imageInVtk->GetScalarComponentAsDouble(testPoint3DInIndex[0], testPoint3DInIndex[1], // testPoint3DInIndex[2], 0); /* Test that everything is working equally if vtkoutput is used instead of the default output * from mitk ImageToImageFilter */ mitk::ExtractSliceFilter::Pointer slicerWithVtkOutput = mitk::ExtractSliceFilter::New(); slicerWithVtkOutput->SetInput(imageInMitk); slicerWithVtkOutput->SetWorldGeometry(plane); slicerWithVtkOutput->SetVtkOutputRequest(true); slicerWithVtkOutput->Update(); vtkSmartPointer vtkImageByVtkOutput = slicerWithVtkOutput->GetVtkOutput(); double PixelvalueByVtkOutput = vtkImageByVtkOutput->GetScalarComponentAsDouble(n1, n2, 0, 0); MITK_TEST_CONDITION(PixelvalueByMitkOutput == PixelvalueByVtkOutput, "testing conversion of image output vtk->mitk by reslicer"); /*================ log outputs ===========================*/ #ifdef EXTRACTOR_DEBUG MITK_INFO << "\n" << "TESTINFO index: " << sliceindex << " orientation: " << orientation << " frontside: " << isFrontside << " rotated: " << isRotated; MITK_INFO << "\n" << "slice index to world: " << sliceIndexToWorld; MITK_INFO << "\n" << "image index to world: " << imageIndexToWorld; MITK_INFO << "\n" << "vtk: slice: " << PixelvalueByMitkOutput << ", image: " << valueVTKinImage; MITK_INFO << "\n" << "testPoint3D InWorld" << testPoint3DInWorld << " is " << testPoint2DInIndex << " in 2D"; MITK_INFO << "\n" << "randoms: " << ((int)n1) << ", " << ((int)n2); MITK_INFO << "\n" << "point is inside plane: " << plane->IsInside(testPoint3DInWorld) << " and volume: " << imageInMitk->GetGeometry()->IsInside(testPoint3DInWorld); MITK_INFO << "\n" << "volume idx: " << testPoint3DInIndex << " = " << valueAt3DVolume; MITK_INFO << "\n" << "volume world: " << testPoint3DInWorld << " = " << valueAt3DVolumeByWorld; MITK_INFO << "\n" << "slice idx: " << testPoint2DInIndex << " = " << valueAtSlice; itk::Index<3> curr; curr[0] = curr[1] = curr[2] = 0; for (int i = 0; i < 32; ++i) { for (int j = 0; j < 32; ++j) { ++curr[1]; if (SliceReadAccessor.GetPixelByIndex(curr) == valueAt3DVolume) { MITK_INFO << "\n" << valueAt3DVolume << " MATCHED mitk " << curr; } } curr[1] = 0; ++curr[0]; } typedef itk::Image Image2DType; Image2DType::Pointer img = Image2DType::New(); CastToItkImage(slice, img); typedef itk::ImageRegionConstIterator Iterator2D; Iterator2D iter(img, img->GetLargestPossibleRegion()); iter.GoToBegin(); while (!iter.IsAtEnd()) { if (img->GetPixel(iter.GetIndex()) == valueAt3DVolume) MITK_INFO << "\n" << valueAt3DVolume << " MATCHED itk " << iter.GetIndex(); ++iter; } #endif // EXTRACTOR_DEBUG } /* random a float value */ static float randFloat() { return (((float)rand() + 1.0) / ((float)RAND_MAX + 1.0)) + (((float)rand() + 1.0) / ((float)RAND_MAX + 1.0)) / ((float)RAND_MAX + 1.0); } /* create a sphere with the size of the given testVolumeSize*/ static void InitializeTestVolume() { #ifdef CREATE_VOLUME // do sphere creation ItkVolumeGeneration(); #ifdef SAVE_VOLUME // save in file mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New(); writer->SetInput(TestVolume); std::string file; std::ostringstream filename; filename << "C:\\home\\schroedt\\MITK\\Modules\\ImageExtraction\\Testing\\Data\\sphere_"; filename << TestvolumeSize; filename << ".nrrd"; file = filename.str(); writer->SetFileName(file); writer->Update(); #endif // SAVE_VOLUME #endif #ifndef CREATE_VOLUME // read from file mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); std::string filename = locator->FindFile("sphere_512.nrrd.mhd", "Modules/ImageExtraction/Testing/Data"); TestVolume = mitk::IOUtil::Load(filename); #endif #ifdef CALC_TESTFAILURE_DEVIATION // get the TestFailureDeviation in % AccessFixedDimensionByItk(TestVolume, CalcTestFailureDeviation, 3); #endif } // the test result of the sphere reslice struct SliceProperties { double planeDistanceToSphereCenter; double diameterInMM; double diameterInPixel; double diameterCalculated; double percentageRadiusToPixel; double areaCalculated; double areaInPixel; double percentageAreaCalcToPixel; }; static mitk::Image::Pointer TestVolume; static double TestvolumeSize; static mitk::PlaneGeometry::Pointer TestPlane; static std::string TestName; static unsigned char pixelValueSet; static SliceProperties testResults; static double TestFailureDeviation; private: /* * Generate a sphere with a radius of TestvolumeSize / 4.0 */ static void ItkVolumeGeneration() { typedef itk::Image TestVolumeType; typedef itk::ImageRegionConstIterator ImageIterator; TestVolumeType::Pointer sphereImage = TestVolumeType::New(); TestVolumeType::IndexType start; start[0] = start[1] = start[2] = 0; TestVolumeType::SizeType size; size[0] = size[1] = size[2] = TestvolumeSize; TestVolumeType::RegionType imgRegion; imgRegion.SetSize(size); imgRegion.SetIndex(start); sphereImage->SetRegions(imgRegion); - sphereImage->SetSpacing(1.0); + sphereImage->SetSpacing(mitk::Vector(1.0)); sphereImage->Allocate(); sphereImage->FillBuffer(0); mitk::Vector3D center; center[0] = center[1] = center[2] = TestvolumeSize / 2.0; double radius = TestvolumeSize / 4.0; double pixelValue = pixelValueSet; ImageIterator imageIterator(sphereImage, sphereImage->GetLargestPossibleRegion()); imageIterator.GoToBegin(); mitk::Vector3D currentVoxelInIndex; while (!imageIterator.IsAtEnd()) { currentVoxelInIndex[0] = imageIterator.GetIndex()[0]; currentVoxelInIndex[1] = imageIterator.GetIndex()[1]; currentVoxelInIndex[2] = imageIterator.GetIndex()[2]; double distanceToCenter = (center + (currentVoxelInIndex * -1.0)).GetNorm(); // if distance to center is smaller then the radius of the sphere if (distanceToCenter < radius) { sphereImage->SetPixel(imageIterator.GetIndex(), pixelValue); } ++imageIterator; } CastToMitkImage(sphereImage, TestVolume); } /* calculate the deviation of the voxel object to the mathematical sphere object. * this is use to make a statement about the accuracy of the resliced image, eg. the circle's diameter or area. */ template static void CalcTestFailureDeviation(itk::Image *inputImage) { typedef itk::Image InputImageType; typedef itk::ImageRegionConstIterator ImageIterator; ImageIterator iterator(inputImage, inputImage->GetLargestPossibleRegion()); iterator.GoToBegin(); int volumeInPixel = 0; while (!iterator.IsAtEnd()) { if (inputImage->GetPixel(iterator.GetIndex()) == pixelValueSet) volumeInPixel++; ++iterator; } double diameter = TestvolumeSize / 2.0; double volumeCalculated = (1.0 / 6.0) * 3.14159265358979 * std::pow(diameter, 3); double volumeDeviation = std::abs(100 - (100 / volumeCalculated * volumeInPixel)); typename InputImageType::IndexType index; index[0] = index[1] = TestvolumeSize / 2.0; index[2] = 0; int sumpixels = 0; while (index[2] < TestvolumeSize) { if (inputImage->GetPixel(index) == pixelValueSet) sumpixels++; index[2] += 1; } double diameterDeviation = std::abs(100 - (100 / diameter * sumpixels)); #ifdef DEBUG MITK_INFO << "volume deviation: " << volumeDeviation << " diameter deviation:" << diameterDeviation; #endif mitkExtractSliceFilterTestClass::TestFailureDeviation = (volumeDeviation + diameterDeviation) / 2.0; } }; /*================ #END class ================*/ /*================#BEGIN Instantiation of members ================*/ mitk::Image::Pointer mitkExtractSliceFilterTestClass::TestVolume = mitk::Image::New(); double mitkExtractSliceFilterTestClass::TestvolumeSize = 256.0; mitk::PlaneGeometry::Pointer mitkExtractSliceFilterTestClass::TestPlane = mitk::PlaneGeometry::New(); std::string mitkExtractSliceFilterTestClass::TestName = ""; unsigned char mitkExtractSliceFilterTestClass::pixelValueSet = 255; mitkExtractSliceFilterTestClass::SliceProperties mitkExtractSliceFilterTestClass::testResults = { -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0}; double mitkExtractSliceFilterTestClass::TestFailureDeviation = 0.0; /*================ #END Instantiation of members ================*/ /*================ #BEGIN test main ================*/ int mitkExtractSliceFilterTest(int /*argc*/, char * /*argv*/ []) { MITK_TEST_BEGIN("mitkExtractSliceFilterTest") // pixelvalue based testing mitkExtractSliceFilterTestClass::PixelvalueBasedTest(); // initialize sphere test volume mitkExtractSliceFilterTestClass::InitializeTestVolume(); mitk::Vector3D spacing = mitkExtractSliceFilterTestClass::TestVolume->GetGeometry()->GetSpacing(); // the center of the sphere = center of image double sphereCenter = mitkExtractSliceFilterTestClass::TestvolumeSize / 2.0; double planeSize = mitkExtractSliceFilterTestClass::TestvolumeSize; /* axial plane */ mitk::PlaneGeometry::Pointer geometryAxial = mitk::PlaneGeometry::New(); geometryAxial->InitializeStandardPlane( planeSize, planeSize, spacing, mitk::AnatomicalPlane::Axial, sphereCenter, false, true); geometryAxial->ChangeImageGeometryConsideringOriginOffset(true); mitk::Point3D origin = geometryAxial->GetOrigin(); mitk::Vector3D normal; normal = geometryAxial->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 // geometryAxial->SetOrigin(origin); mitkExtractSliceFilterTestClass::TestSlice(geometryAxial, "Testing axial plane"); /* end axial plane */ /* sagittal plane */ mitk::PlaneGeometry::Pointer geometrySagittal = mitk::PlaneGeometry::New(); geometrySagittal->InitializeStandardPlane( planeSize, planeSize, spacing, mitk::AnatomicalPlane::Sagittal, sphereCenter, true, false); geometrySagittal->ChangeImageGeometryConsideringOriginOffset(true); origin = geometrySagittal->GetOrigin(); normal = geometrySagittal->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 // geometrySagittal->SetOrigin(origin); mitkExtractSliceFilterTestClass::TestSlice(geometrySagittal, "Testing sagittal plane"); /* sagittal plane */ /* sagittal shifted plane */ mitk::PlaneGeometry::Pointer geometrySagittalShifted = mitk::PlaneGeometry::New(); geometrySagittalShifted->InitializeStandardPlane( planeSize, planeSize, spacing, mitk::AnatomicalPlane::Sagittal, (sphereCenter - 14), true, false); geometrySagittalShifted->ChangeImageGeometryConsideringOriginOffset(true); origin = geometrySagittalShifted->GetOrigin(); normal = geometrySagittalShifted->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 // geometrySagittalShifted->SetOrigin(origin); mitkExtractSliceFilterTestClass::TestSlice(geometrySagittalShifted, "Testing sagittal plane shifted"); /* end sagittal shifted plane */ /* coronal plane */ mitk::PlaneGeometry::Pointer geometryCoronal = mitk::PlaneGeometry::New(); geometryCoronal->InitializeStandardPlane( planeSize, planeSize, spacing, mitk::AnatomicalPlane::Coronal, sphereCenter, true, false); geometryCoronal->ChangeImageGeometryConsideringOriginOffset(true); origin = geometryCoronal->GetOrigin(); normal = geometryCoronal->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 // geometryCoronal->SetOrigin(origin); mitkExtractSliceFilterTestClass::TestSlice(geometryCoronal, "Testing coronal plane"); /* end coronal plane */ /* oblique plane */ mitk::PlaneGeometry::Pointer obliquePlane = mitk::PlaneGeometry::New(); obliquePlane->InitializeStandardPlane( planeSize, planeSize, spacing, mitk::AnatomicalPlane::Sagittal, sphereCenter, true, false); obliquePlane->ChangeImageGeometryConsideringOriginOffset(true); origin = obliquePlane->GetOrigin(); normal = obliquePlane->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 // obliquePlane->SetOrigin(origin); mitk::Vector3D rotationVector; rotationVector[0] = 0.2; rotationVector[1] = 0.4; rotationVector[2] = 0.62; float degree = 37.0; mitk::RotationOperation *op = new mitk::RotationOperation(mitk::OpROTATE, obliquePlane->GetCenter(), rotationVector, degree); obliquePlane->ExecuteOperation(op); delete op; mitkExtractSliceFilterTestClass::TestSlice(obliquePlane, "Testing oblique plane"); /* end oblique plane */ #ifdef SHOW_SLICE_IN_RENDER_WINDOW /*================ #BEGIN vtk render code ================*/ // set reslicer for renderwindow mitk::Image::Pointer pic = mitk::IOUtil::Load(filename); vtkSmartPointer slicer = vtkSmartPointer::New(); slicer->SetInput(pic->GetVtkImageData()); mitk::PlaneGeometry::Pointer obliquePl = mitk::PlaneGeometry::New(); obliquePl->InitializeStandardPlane( pic->GetGeometry(), mitk::AnatomicalPlane::Sagittal, pic->GetGeometry()->GetCenter()[0], true, false); obliquePl->ChangeImageGeometryConsideringOriginOffset(true); mitk::Point3D origin2 = obliquePl->GetOrigin(); mitk::Vector3D n; n = obliquePl->GetNormal(); n.Normalize(); origin2 += n * 0.5; // pixelspacing is 1, so half the spacing is 0.5 obliquePl->SetOrigin(origin2); mitk::Vector3D rotation; rotation[0] = 0.534307; rotation[1] = 0.000439605; rotation[2] = 0.423017; MITK_INFO << rotation; mitk::RotationOperation *operation = new mitk::RotationOperation(mitk::OpROTATE, obliquePl->GetCenter(), rotationVector, degree); obliquePl->ExecuteOperation(operation); delete operation; double origin[3]; origin[0] = obliquePl->GetOrigin()[0]; origin[1] = obliquePl->GetOrigin()[1]; origin[2] = obliquePl->GetOrigin()[2]; slicer->SetResliceAxesOrigin(origin); mitk::Vector3D right, bottom, normal; right = obliquePl->GetAxisVector(0); bottom = obliquePl->GetAxisVector(1); normal = obliquePl->GetNormal(); right.Normalize(); bottom.Normalize(); normal.Normalize(); double cosines[9]; mitk::vnl2vtk(right.GetVnlVector(), cosines); // x mitk::vnl2vtk(bottom.GetVnlVector(), cosines + 3); // y mitk::vnl2vtk(normal.GetVnlVector(), cosines + 6); // n slicer->SetResliceAxesDirectionCosines(cosines); slicer->SetOutputDimensionality(2); slicer->Update(); // set vtk renderwindow vtkSmartPointer vtkPlane = vtkSmartPointer::New(); vtkPlane->SetOrigin(0.0, 0.0, 0.0); // These two points define the axes of the plane in combination with the origin. // Point 1 is the x-axis and point 2 the y-axis. // Each plane is transformed according to the view (axial, coronal and sagittal) afterwards. vtkPlane->SetPoint1(1.0, 0.0, 0.0); // P1: (xMax, yMin, depth) vtkPlane->SetPoint2(0.0, 1.0, 0.0); // P2: (xMin, yMax, depth) // these are not the correct values for all slices, only a square plane by now vtkSmartPointer imageMapper = vtkSmartPointer::New(); imageMapper->SetInputConnection(vtkPlane->GetOutputPort()); vtkSmartPointer lookupTable = vtkSmartPointer::New(); // built a default lookuptable lookupTable->SetRampToLinear(); lookupTable->SetSaturationRange(0.0, 0.0); lookupTable->SetHueRange(0.0, 0.0); lookupTable->SetValueRange(0.0, 1.0); lookupTable->Build(); // map all black values to transparent lookupTable->SetTableValue(0, 0.0, 0.0, 0.0, 0.0); lookupTable->SetRange(-255.0, 255.0); // lookupTable->SetRange(-1022.0, 1184.0);//pic3D range vtkSmartPointer texture = vtkSmartPointer::New(); texture->SetInput(slicer->GetOutput()); texture->SetLookupTable(lookupTable); texture->SetMapColorScalarsThroughLookupTable(true); vtkSmartPointer imageActor = vtkSmartPointer::New(); imageActor->SetMapper(imageMapper); imageActor->SetTexture(texture); // Setup renderers vtkSmartPointer renderer = vtkSmartPointer::New(); renderer->AddActor(imageActor); // Setup render window vtkSmartPointer renderWindow = vtkSmartPointer::New(); renderWindow->AddRenderer(renderer); // Setup render window interactor vtkSmartPointer renderWindowInteractor = vtkSmartPointer::New(); vtkSmartPointer style = vtkSmartPointer::New(); renderWindowInteractor->SetInteractorStyle(style); // Render and start interaction renderWindowInteractor->SetRenderWindow(renderWindow); // renderer->AddViewProp(imageActor); renderWindow->Render(); renderWindowInteractor->Start(); // always end with this! /*================ #END vtk render code ================*/ #endif // SHOW_SLICE_IN_RENDER_WINDOW MITK_TEST_END() } diff --git a/Modules/Core/test/mitkGeometry3DTest.cpp b/Modules/Core/test/mitkGeometry3DTest.cpp index caa5fc4d1a..c4395e1a28 100644 --- a/Modules/Core/test/mitkGeometry3DTest.cpp +++ b/Modules/Core/test/mitkGeometry3DTest.cpp @@ -1,622 +1,622 @@ /*============================================================================ 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 "mitkGeometry3D.h" #include #include #include "mitkInteractionConst.h" #include "mitkRotationOperation.h" #include #include #include "mitkTestingMacros.h" #include #include bool testGetAxisVectorVariants(mitk::Geometry3D *geometry) { int direction; for (direction = 0; direction < 3; ++direction) { mitk::Vector3D frontToBack(0.); switch (direction) { case 0: frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(true, false, false); break; // 7-3 case 1: frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(false, true, false); break; // 7-5 case 2: frontToBack = geometry->GetCornerPoint(false, false, false) - geometry->GetCornerPoint(false, false, true); break; // 7-2 } std::cout << "Testing GetAxisVector(int) vs GetAxisVector(bool, bool, bool): "; if (mitk::Equal(geometry->GetAxisVector(direction), frontToBack) == false) { std::cout << "[FAILED]" << std::endl; return false; } std::cout << "[PASSED]" << std::endl; } return true; } bool testGetAxisVectorExtent(mitk::Geometry3D *geometry) { int direction; for (direction = 0; direction < 3; ++direction) { if (mitk::Equal(geometry->GetAxisVector(direction).GetNorm(), geometry->GetExtentInMM(direction)) == false) { std::cout << "[FAILED]" << std::endl; return false; } std::cout << "[PASSED]" << std::endl; } return true; } // a part of the test requires axis-parallel coordinates int testIndexAndWorldConsistency(mitk::Geometry3D *geometry3d) { MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems: "); mitk::Point3D origin = geometry3d->GetOrigin(); mitk::Point3D dummy; MITK_TEST_OUTPUT(<< " Testing index->world->index conversion consistency"); geometry3d->WorldToIndex(origin, dummy); geometry3d->IndexToWorld(dummy, dummy); MITK_TEST_CONDITION_REQUIRED(dummy == origin, ""); MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin, mitk::Point3D)==(0,0,0)"); mitk::Point3D globalOrigin; mitk::FillVector3D(globalOrigin, 0, 0, 0); mitk::Point3D originContinuousIndex; geometry3d->WorldToIndex(origin, originContinuousIndex); MITK_TEST_CONDITION_REQUIRED(originContinuousIndex == globalOrigin, ""); MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin, itk::Index)==(0,0,0)"); itk::Index<3> itkindex; geometry3d->WorldToIndex(origin, itkindex); itk::Index<3> globalOriginIndex; mitk::vtk2itk(globalOrigin, globalOriginIndex); MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin-0.5*spacing, itk::Index)==(0,0,0)"); mitk::Vector3D halfSpacingStep = geometry3d->GetSpacing() * 0.5; mitk::Matrix3D rotation; mitk::Point3D originOffCenter = origin - halfSpacingStep; geometry3d->WorldToIndex(originOffCenter, itkindex); MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0)"); originOffCenter = origin + halfSpacingStep; - originOffCenter -= 0.0001; + originOffCenter -= mitk::Vector(0.0001); geometry3d->WorldToIndex(originOffCenter, itkindex); MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); MITK_TEST_OUTPUT(<< " Testing WorldToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)"); originOffCenter = origin + halfSpacingStep; itk::Index<3> global111; mitk::FillVector3D(global111, 1, 1, 1); geometry3d->WorldToIndex(originOffCenter, itkindex); MITK_TEST_CONDITION_REQUIRED(itkindex == global111, ""); MITK_TEST_OUTPUT(<< " Testing WorldToIndex(GetCenter())==BoundingBox.GetCenter: "); mitk::Point3D center = geometry3d->GetCenter(); mitk::Point3D centerContIndex; geometry3d->WorldToIndex(center, centerContIndex); mitk::BoundingBox::ConstPointer boundingBox = geometry3d->GetBoundingBox(); mitk::BoundingBox::PointType centerBounds = boundingBox->GetCenter(); // itk assumes corner based geometry. If our geometry is center based (imageGoe == true), everything needs to be // shifted if (geometry3d->GetImageGeometry()) { centerBounds[0] -= 0.5; centerBounds[1] -= 0.5; centerBounds[2] -= 0.5; } MITK_TEST_CONDITION_REQUIRED(mitk::Equal(centerContIndex, centerBounds), ""); MITK_TEST_OUTPUT(<< " Testing GetCenter()==IndexToWorld(BoundingBox.GetCenter): "); center = geometry3d->GetCenter(); mitk::Point3D centerBoundsInWorldCoords; geometry3d->IndexToWorld(centerBounds, centerBoundsInWorldCoords); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(center, centerBoundsInWorldCoords), ""); return EXIT_SUCCESS; } int testIndexAndWorldConsistencyForVectors(mitk::Geometry3D *geometry3d) { MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems for vectors: "); mitk::Vector3D xAxisMM = geometry3d->GetAxisVector(0); mitk::Vector3D xAxisContinuousIndex; mitk::Vector3D xAxisContinuousIndexDeprecated; mitk::Point3D p, pIndex, origin; origin = geometry3d->GetOrigin(); p[0] = xAxisMM[0]; p[1] = xAxisMM[1]; p[2] = xAxisMM[2]; geometry3d->WorldToIndex(p, pIndex); geometry3d->WorldToIndex(xAxisMM, xAxisContinuousIndexDeprecated); geometry3d->WorldToIndex(xAxisMM, xAxisContinuousIndex); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == pIndex[0], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == pIndex[1], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == pIndex[2], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == xAxisContinuousIndexDeprecated[0], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == xAxisContinuousIndexDeprecated[1], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == xAxisContinuousIndexDeprecated[2], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[0] == pIndex[0], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[1] == pIndex[1], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[2] == pIndex[2], ""); geometry3d->IndexToWorld(xAxisContinuousIndex, xAxisContinuousIndex); geometry3d->IndexToWorld(xAxisContinuousIndexDeprecated, xAxisContinuousIndexDeprecated); geometry3d->IndexToWorld(pIndex, p); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex == xAxisMM, ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[0] == p[0], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[1] == p[1], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndex[2] == p[2], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated == xAxisMM, ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[0] == p[0], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[1] == p[1], ""); MITK_TEST_CONDITION_REQUIRED(xAxisContinuousIndexDeprecated[2] == p[2], ""); return EXIT_SUCCESS; } int testIndexAndWorldConsistencyForIndex(mitk::Geometry3D *geometry3d) { MITK_TEST_OUTPUT(<< "Testing consistency of index and world coordinate systems: "); // creating testing data itk::Index<4> itkIndex4, itkIndex4b; itk::Index<3> itkIndex3, itkIndex3b; itk::Index<2> itkIndex2, itkIndex2b; itk::Index<3> mitkIndex, mitkIndexb; itkIndex4[0] = itkIndex4[1] = itkIndex4[2] = itkIndex4[3] = 4; itkIndex3[0] = itkIndex3[1] = itkIndex3[2] = 6; itkIndex2[0] = itkIndex2[1] = 2; mitkIndex[0] = mitkIndex[1] = mitkIndex[2] = 13; // check for consistency mitk::Point3D point; geometry3d->IndexToWorld(itkIndex2, point); geometry3d->WorldToIndex(point, itkIndex2b); MITK_TEST_CONDITION_REQUIRED(((itkIndex2b[0] == itkIndex2[0]) && (itkIndex2b[1] == itkIndex2[1])), "Testing itk::index<2> for IndexToWorld/WorldToIndex consistency"); geometry3d->IndexToWorld(itkIndex3, point); geometry3d->WorldToIndex(point, itkIndex3b); MITK_TEST_CONDITION_REQUIRED( ((itkIndex3b[0] == itkIndex3[0]) && (itkIndex3b[1] == itkIndex3[1]) && (itkIndex3b[2] == itkIndex3[2])), "Testing itk::index<3> for IndexToWorld/WorldToIndex consistency"); geometry3d->IndexToWorld(itkIndex4, point); geometry3d->WorldToIndex(point, itkIndex4b); MITK_TEST_CONDITION_REQUIRED(((itkIndex4b[0] == itkIndex4[0]) && (itkIndex4b[1] == itkIndex4[1]) && (itkIndex4b[2] == itkIndex4[2]) && (itkIndex4b[3] == 0)), "Testing itk::index<3> for IndexToWorld/WorldToIndex consistency"); geometry3d->IndexToWorld(mitkIndex, point); geometry3d->WorldToIndex(point, mitkIndexb); MITK_TEST_CONDITION_REQUIRED( ((mitkIndexb[0] == mitkIndex[0]) && (mitkIndexb[1] == mitkIndex[1]) && (mitkIndexb[2] == mitkIndex[2])), "Testing mitk::Index for IndexToWorld/WorldToIndex consistency"); return EXIT_SUCCESS; } #include int testItkImageIsCenterBased() { MITK_TEST_OUTPUT(<< "Testing whether itk::Image coordinates are center-based."); typedef itk::Image ItkIntImage3D; ItkIntImage3D::Pointer itkintimage = ItkIntImage3D::New(); ItkIntImage3D::SizeType size; size.Fill(10); mitk::Point3D origin; mitk::FillVector3D(origin, 2, 3, 7); itkintimage->Initialize(); itkintimage->SetRegions(size); itkintimage->SetOrigin(origin); std::cout << "[PASSED]" << std::endl; MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToContinuousIndex(origin)==(0,0,0)"); mitk::Point3D globalOrigin; mitk::FillVector3D(globalOrigin, 0, 0, 0); itk::ContinuousIndex originContinuousIndex; itkintimage->TransformPhysicalPointToContinuousIndex(origin, originContinuousIndex); MITK_TEST_CONDITION_REQUIRED(originContinuousIndex == globalOrigin, ""); MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToIndex(origin)==(0,0,0)"); itk::Index<3> itkindex; itkintimage->TransformPhysicalPointToIndex(origin, itkindex); itk::Index<3> globalOriginIndex; mitk::vtk2itk(globalOrigin, globalOriginIndex); MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToIndex(origin-0.5*spacing)==(0,0,0)"); mitk::Vector3D halfSpacingStep = itkintimage->GetSpacing() * 0.5; mitk::Matrix3D rotation; mitk::Point3D originOffCenter = origin - halfSpacingStep; itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex); MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); MITK_TEST_OUTPUT( << " Testing itk::Image::TransformPhysicalPointToIndex(origin+0.5*spacing-eps, itk::Index)==(0,0,0)"); originOffCenter = origin + halfSpacingStep; - originOffCenter -= 0.0001; + originOffCenter -= mitk::Vector(0.0001); itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex); MITK_TEST_CONDITION_REQUIRED(itkindex == globalOriginIndex, ""); MITK_TEST_OUTPUT(<< " Testing itk::Image::TransformPhysicalPointToIndex(origin+0.5*spacing, itk::Index)==(1,1,1)"); originOffCenter = origin + halfSpacingStep; itk::Index<3> global111; mitk::FillVector3D(global111, 1, 1, 1); itkintimage->TransformPhysicalPointToIndex(originOffCenter, itkindex); MITK_TEST_CONDITION_REQUIRED(itkindex == global111, ""); MITK_TEST_OUTPUT(<< "=> Yes, itk::Image coordinates are center-based."); return EXIT_SUCCESS; } int testGeometry3D(bool imageGeometry) { // Build up a new image Geometry mitk::Geometry3D::Pointer geometry3d = mitk::Geometry3D::New(); float bounds[] = {-10.0, 17.0, -12.0, 188.0, 13.0, 211.0}; MITK_TEST_OUTPUT(<< "Initializing"); geometry3d->Initialize(); MITK_TEST_OUTPUT(<< "Setting ImageGeometry to " << imageGeometry); geometry3d->SetImageGeometry(imageGeometry); MITK_TEST_OUTPUT(<< "Setting bounds by SetFloatBounds(): " << bounds); geometry3d->SetFloatBounds(bounds); MITK_TEST_OUTPUT(<< "Testing AxisVectors"); if (testGetAxisVectorVariants(geometry3d) == false) return EXIT_FAILURE; if (testGetAxisVectorExtent(geometry3d) == false) return EXIT_FAILURE; MITK_TEST_OUTPUT(<< "Creating an AffineTransform3D transform"); mitk::AffineTransform3D::MatrixType matrix; matrix.SetIdentity(); matrix(1, 1) = 2; mitk::AffineTransform3D::Pointer transform; transform = mitk::AffineTransform3D::New(); transform->SetMatrix(matrix); MITK_TEST_OUTPUT(<< "Testing a SetIndexToWorldTransform"); geometry3d->SetIndexToWorldTransform(transform); MITK_TEST_OUTPUT(<< "Testing correctness of value returned by GetSpacing"); const mitk::Vector3D &spacing1 = geometry3d->GetSpacing(); mitk::Vector3D expectedSpacing; expectedSpacing.Fill(1.0); expectedSpacing[1] = 2; if (mitk::Equal(spacing1, expectedSpacing) == false) { MITK_TEST_OUTPUT(<< " [FAILED]"); return EXIT_FAILURE; } MITK_TEST_OUTPUT(<< "Testing a Compose(transform)"); geometry3d->Compose(transform); MITK_TEST_OUTPUT(<< "Testing correctness of value returned by GetSpacing"); const mitk::Vector3D &spacing2 = geometry3d->GetSpacing(); expectedSpacing[1] = 4; if (mitk::Equal(spacing2, expectedSpacing) == false) { MITK_TEST_OUTPUT(<< " [FAILED]"); return EXIT_FAILURE; } MITK_TEST_OUTPUT(<< "Testing correctness of SetSpacing"); mitk::Vector3D newspacing; mitk::FillVector3D(newspacing, 1.5, 2.5, 3.5); geometry3d->SetSpacing(newspacing); const mitk::Vector3D &spacing3 = geometry3d->GetSpacing(); if (mitk::Equal(spacing3, newspacing) == false) { MITK_TEST_OUTPUT(<< " [FAILED]"); return EXIT_FAILURE; } // Separate Test function for Index and World consistency testIndexAndWorldConsistency(geometry3d); testIndexAndWorldConsistencyForVectors(geometry3d); testIndexAndWorldConsistencyForIndex(geometry3d); MITK_TEST_OUTPUT(<< "Testing a rotation of the geometry"); double angle = 35.0; mitk::Vector3D rotationVector; mitk::FillVector3D(rotationVector, 1, 0, 0); mitk::Point3D center = geometry3d->GetCenter(); auto op = new mitk::RotationOperation(mitk::OpROTATE, center, rotationVector, angle); geometry3d->ExecuteOperation(op); MITK_TEST_OUTPUT(<< "Testing mitk::GetRotation() and success of rotation"); mitk::Matrix3D rotation; mitk::GetRotation(geometry3d, rotation); mitk::Vector3D voxelStep = rotation * newspacing; mitk::Vector3D voxelStepIndex; geometry3d->WorldToIndex(voxelStep, voxelStepIndex); mitk::Vector3D expectedVoxelStepIndex; expectedVoxelStepIndex.Fill(1); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(voxelStepIndex, expectedVoxelStepIndex), ""); delete op; std::cout << "[PASSED]" << std::endl; MITK_TEST_OUTPUT(<< "Testing that ImageGeometry is still " << imageGeometry); MITK_TEST_CONDITION_REQUIRED(geometry3d->GetImageGeometry() == imageGeometry, ""); // Test if the translate function moves the origin correctly. mitk::Point3D oldOrigin = geometry3d->GetOrigin(); // use some random values for translation mitk::Vector3D translationVector; translationVector.SetElement(0, 17.5f); translationVector.SetElement(1, -32.3f); translationVector.SetElement(2, 4.0f); // compute ground truth mitk::Point3D tmpResult = geometry3d->GetOrigin() + translationVector; geometry3d->Translate(translationVector); MITK_TEST_CONDITION(mitk::Equal(geometry3d->GetOrigin(), tmpResult), "Testing if origin was translated."); translationVector *= -1; // vice versa geometry3d->Translate(translationVector); MITK_TEST_CONDITION(mitk::Equal(geometry3d->GetOrigin(), oldOrigin), "Testing if the translation could be done vice versa."); return EXIT_SUCCESS; } int testGeometryAfterCasting() { // Epsilon. Allowed difference for rotationvalue float eps = 0.0001; // Cast ITK and MITK images and see if geometry stays typedef itk::Image Image2DType; typedef itk::Image Image3DType; // Create 3D ITK Image from Scratch, cast to 3D MITK image, compare Geometries Image3DType::Pointer image3DItk = Image3DType::New(); Image3DType::RegionType myRegion; Image3DType::SizeType mySize; Image3DType::IndexType myIndex; Image3DType::SpacingType mySpacing; Image3DType::DirectionType myDirection, rotMatrixX, rotMatrixY, rotMatrixZ; mySpacing[0] = 31; mySpacing[1] = 0.1; mySpacing[2] = 2.9; myIndex[0] = -15; myIndex[1] = 15; myIndex[2] = 12; mySize[0] = 10; mySize[1] = 2; mySize[2] = 5; myRegion.SetSize(mySize); myRegion.SetIndex(myIndex); image3DItk->SetSpacing(mySpacing); image3DItk->SetRegions(myRegion); image3DItk->Allocate(); image3DItk->FillBuffer(0); myDirection.SetIdentity(); rotMatrixX.SetIdentity(); rotMatrixY.SetIdentity(); rotMatrixZ.SetIdentity(); mitk::Image::Pointer mitkImage; // direction [row] [column] MITK_TEST_OUTPUT(<< "Casting a rotated 3D ITK Image to a MITK Image and check if Geometry is still same"); for (double rotX = 0; rotX < itk::Math::pi * 2; rotX += itk::Math::pi * 0.4) { // Set Rotation X rotMatrixX[1][1] = cos(rotX); rotMatrixX[1][2] = -sin(rotX); rotMatrixX[2][1] = sin(rotX); rotMatrixX[2][2] = cos(rotX); for (double rotY = 0; rotY < itk::Math::pi * 2; rotY += itk::Math::pi * 0.3) { // Set Rotation Y rotMatrixY[0][0] = cos(rotY); rotMatrixY[0][2] = sin(rotY); rotMatrixY[2][0] = -sin(rotY); rotMatrixY[2][2] = cos(rotY); for (double rotZ = 0; rotZ < itk::Math::pi * 2; rotZ += itk::Math::pi * 0.2) { // Set Rotation Z rotMatrixZ[0][0] = cos(rotZ); rotMatrixZ[0][1] = -sin(rotZ); rotMatrixZ[1][0] = sin(rotZ); rotMatrixZ[1][1] = cos(rotZ); // Multiply matrizes myDirection = myDirection * rotMatrixX * rotMatrixY * rotMatrixZ; image3DItk->SetDirection(myDirection); mitk::CastToMitkImage(image3DItk, mitkImage); const mitk::AffineTransform3D::MatrixType &matrix = mitkImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix(); for (int row = 0; row < 3; row++) { for (int col = 0; col < 3; col++) { double mitkValue = matrix[row][col] / mitkImage->GetGeometry()->GetSpacing()[col]; double itkValue = myDirection[row][col]; double diff = mitkValue - itkValue; // if you decrease this value, you can see that there might be QUITE high inaccuracy!!! if (diff > eps) // need to check, how exact it SHOULD be .. since it is NOT EXACT! { std::cout << "Had a difference of : " << diff; std::cout << "Error: Casting altered Geometry!"; std::cout << "ITK Matrix:\n" << myDirection; std::cout << "Mitk Matrix (With Spacing):\n" << matrix; std::cout << "Mitk Spacing: " << mitkImage->GetGeometry()->GetSpacing(); MITK_TEST_CONDITION_REQUIRED(false == true, ""); return false; } } } } } } // Create 2D ITK Image from Scratch, cast to 2D MITK image, compare Geometries Image2DType::Pointer image2DItk = Image2DType::New(); Image2DType::RegionType myRegion2D; Image2DType::SizeType mySize2D; Image2DType::IndexType myIndex2D; Image2DType::SpacingType mySpacing2D; Image2DType::DirectionType myDirection2D, rotMatrix; mySpacing2D[0] = 31; mySpacing2D[1] = 0.1; myIndex2D[0] = -15; myIndex2D[1] = 15; mySize2D[0] = 10; mySize2D[1] = 2; myRegion2D.SetSize(mySize2D); myRegion2D.SetIndex(myIndex2D); image2DItk->SetSpacing(mySpacing2D); image2DItk->SetRegions(myRegion2D); image2DItk->Allocate(); image2DItk->FillBuffer(0); myDirection2D.SetIdentity(); rotMatrix.SetIdentity(); // direction [row] [column] MITK_TEST_OUTPUT(<< "Casting a rotated 2D ITK Image to a MITK Image and check if Geometry is still same"); for (double rotTheta = 0; rotTheta < itk::Math::pi * 2; rotTheta += itk::Math::pi * 0.2) { // Set Rotation rotMatrix[0][0] = cos(rotTheta); rotMatrix[0][1] = -sin(rotTheta); rotMatrix[1][0] = sin(rotTheta); rotMatrix[1][1] = cos(rotTheta); // Multiply matrizes myDirection2D = myDirection2D * rotMatrix; image2DItk->SetDirection(myDirection2D); mitk::CastToMitkImage(image2DItk, mitkImage); const mitk::AffineTransform3D::MatrixType &matrix = mitkImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix(); // Compare MITK and ITK matrix for (int row = 0; row < 3; row++) { for (int col = 0; col < 3; col++) { double mitkValue = matrix[row][col] / mitkImage->GetGeometry()->GetSpacing()[col]; if ((row == 2) && (col == row)) { if (mitkValue != 1) { MITK_TEST_OUTPUT(<< "After casting a 2D ITK to 3D MITK images, MITK matrix values for 0|2, 1|2, 2|0, 2|1 " "MUST be 0 and value for 2|2 must be 1"); return false; } } else if ((row == 2) || (col == 2)) { if (mitkValue != 0) { MITK_TEST_OUTPUT(<< "After casting a 2D ITK to 3D MITK images, MITK matrix values for 0|2, 1|2, 2|0, 2|1 " "MUST be 0 and value for 2|2 must be 1"); return false; } } else { double itkValue = myDirection2D[row][col]; double diff = mitkValue - itkValue; // if you decrease this value, you can see that there might be QUITE high inaccuracy!!! if (diff > eps) // need to check, how exact it SHOULD be .. since it is NOT EXACT! { std::cout << "Had a difference of : " << diff; std::cout << "Error: Casting altered Geometry!"; std::cout << "ITK Matrix:\n" << myDirection2D; std::cout << "Mitk Matrix (With Spacing):\n" << matrix; std::cout << "Mitk Spacing: " << mitkImage->GetGeometry()->GetSpacing(); MITK_TEST_CONDITION_REQUIRED(false == true, ""); return false; } } } } } // THIS WAS TESTED: // 2D ITK -> 2D MITK, // 3D ITK -> 3D MITK, // Still need to test: 2D MITK Image with ADDITIONAL INFORMATION IN MATRIX -> 2D ITK // 1. Possibility: 3x3 MITK matrix can be converted without loss into 2x2 ITK matrix // 2. Possibility: 3x3 MITK matrix can only be converted with loss into 2x2 ITK matrix // .. before implementing this, we wait for further development in geometry classes (e.g. Geoemtry3D::SetRotation(..)) return EXIT_SUCCESS; } int mitkGeometry3DTest(int /*argc*/, char * /*argv*/ []) { MITK_TEST_BEGIN(mitkGeometry3DTest); int result; MITK_TEST_CONDITION_REQUIRED((result = testItkImageIsCenterBased()) == EXIT_SUCCESS, ""); MITK_TEST_OUTPUT(<< "Running main part of test with ImageGeometry = false"); MITK_TEST_CONDITION_REQUIRED((result = testGeometry3D(false)) == EXIT_SUCCESS, ""); MITK_TEST_OUTPUT(<< "Running main part of test with ImageGeometry = true"); MITK_TEST_CONDITION_REQUIRED((result = testGeometry3D(true)) == EXIT_SUCCESS, ""); MITK_TEST_OUTPUT(<< "Running test to see if Casting MITK to ITK and the other way around destroys geometry"); MITK_TEST_CONDITION_REQUIRED((result = testGeometryAfterCasting()) == EXIT_SUCCESS, ""); MITK_TEST_END(); return EXIT_SUCCESS; } diff --git a/Modules/Core/test/mitkPointTypeConversionTest.cpp b/Modules/Core/test/mitkPointTypeConversionTest.cpp index 9b706c1e03..22d6267c35 100644 --- a/Modules/Core/test/mitkPointTypeConversionTest.cpp +++ b/Modules/Core/test/mitkPointTypeConversionTest.cpp @@ -1,139 +1,139 @@ /*============================================================================ 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 "mitkTestFixture.h" #include "itkPoint.h" #include "mitkNumericConstants.h" #include "mitkNumericTypes.h" // for Equal method #include "mitkPoint.h" #include "mitkTestingMacros.h" #include "vtkPoints.h" #include "vtkSmartPointer.h" #include using namespace mitk; class mitkPointTypeConversionTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPointTypeConversionTestSuite); MITK_TEST(Vector2Point); MITK_TEST(Mitk2Itk_PointCompatibility); MITK_TEST(Itk2Mitk_PointCompatibility); MITK_TEST(Vtk2Mitk_PointCompatibility); MITK_TEST(Mitk2Pod_PointCompatibility); MITK_TEST(Pod2Mitk_PointCompatibility); CPPUNIT_TEST_SUITE_END(); private: vtkSmartPointer a_vtkPoints; ScalarType originalValues[3]; ScalarType valuesToCopy[3]; /** * @brief Convenience method to test if one vector has been assigned successfully to the other. * * More specifically, tests if v1 = v2 was performed correctly. * * @param v1 The vector v1 of the assignment v1 = v2 * @param v2 The vector v2 of the assignment v1 = v2 * @param v1Name The type name of v1 (e.g.: mitk::Vector3D). Necessary for the correct test output. * @param v2Name The type name of v2 (e.g.: mitk::Vector3D). Necessary for the correct test output. * @param eps defines the allowed tolerance when testing for equality. */ template void TestForEquality(T1 v1, T2 v2, std::string v1Name, std::string v2Name, ScalarType eps = mitk::eps) { CPPUNIT_ASSERT_EQUAL_MESSAGE( "\nAssigning " + v2Name + " to " + v1Name + ":\n both are equal", true, EqualArray(v1, v2, 3, eps)); } public: void setUp(void) override { FillVector3D(originalValues, 1.0, 2.0, 3.0); FillVector3D(valuesToCopy, 4.0, 5.0, 6.0); a_vtkPoints = vtkSmartPointer::New(); a_vtkPoints->Initialize(); } void tearDown(void) override { // a_vtkPoints = nullptr; } void Mitk2Itk_PointCompatibility() { mitk::Point3D point3D = valuesToCopy; itk::Point itkPoint3D = point3D; TestForEquality(itkPoint3D, point3D, "itk::Point", "mitk:Point"); } void Itk2Mitk_PointCompatibility() { itk::Point itkPoint3D = valuesToCopy; mitk::Point3D point3D = itkPoint3D; TestForEquality(point3D, itkPoint3D, "mitk:Point", "itk::Point"); } void Vtk2Mitk_PointCompatibility() { a_vtkPoints->InsertNextPoint(valuesToCopy); double vtkPoint[3]; a_vtkPoints->GetPoint(0, vtkPoint); mitk::Point3D point3D = vtkPoint; TestForEquality(point3D, vtkPoint, "mitk:Point", "vtkPoint"); } void Mitk2Pod_PointCompatibility() { ScalarType podPoint[] = {1.0, 2.0, 3.0}; mitk::Point3D point3D = valuesToCopy; point3D.ToArray(podPoint); TestForEquality(podPoint, point3D, "POD point", "mitk::Point"); } void Pod2Mitk_PointCompatibility() { ScalarType podPoint[] = {4.0, 5.0, 6.0}; itk::Point point3D = podPoint; TestForEquality(point3D, podPoint, "mitk::Point3D", "POD point"); } void Vector2Point() { itk::Vector vector3D = originalValues; - itk::Point point3D = vector3D; + itk::Point point3D(vector3D.GetDataPointer()); TestForEquality(point3D, vector3D, "mitk::Point", "mitk::Vector"); } }; MITK_TEST_SUITE_REGISTRATION(mitkPointTypeConversion) diff --git a/Modules/IGT/TrackingDevices/mitkVirtualTrackingDevice.cpp b/Modules/IGT/TrackingDevices/mitkVirtualTrackingDevice.cpp index cdf8a4a56f..37c67f5fbd 100644 --- a/Modules/IGT/TrackingDevices/mitkVirtualTrackingDevice.cpp +++ b/Modules/IGT/TrackingDevices/mitkVirtualTrackingDevice.cpp @@ -1,307 +1,307 @@ /*============================================================================ 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 "mitkVirtualTrackingDevice.h" #include "mitkIGTTimeStamp.h" #include "mitkIGTException.h" #include #include #include #include #include #include mitk::VirtualTrackingDevice::VirtualTrackingDevice() : mitk::TrackingDevice(), m_AllTools(), m_RefreshRate(100), m_NumberOfControlPoints(20), m_GaussianNoiseEnabled(false), m_MeanDistributionParam(0.0), m_DeviationDistributionParam(1.0) { m_Data = mitk::VirtualTrackerTypeInformation::GetDeviceDataVirtualTracker(); m_Bounds[0] = m_Bounds[2] = m_Bounds[4] = -400.0; // initialize bounds to -400 ... +400 (mm) cube m_Bounds[1] = m_Bounds[3] = m_Bounds[5] = 400.0; } mitk::VirtualTrackingDevice::~VirtualTrackingDevice() { if (this->GetState() == Tracking) { this->StopTracking(); } if (this->GetState() == Ready) { this->CloseConnection(); } /* cleanup tracking thread */ if (m_Thread.joinable()) m_Thread.join(); m_AllTools.clear(); } mitk::TrackingTool* mitk::VirtualTrackingDevice::AddTool(const char* toolName) { //if (this->GetState() == Tracking) //{ // return nullptr; //} mitk::VirtualTrackingTool::Pointer t = mitk::VirtualTrackingTool::New(); t->SetToolName(toolName); t->SetVelocity(0.1); this->InitializeSpline(t); std::lock_guard lock(m_ToolsMutex); // lock and unlock the mutex m_AllTools.push_back(t); return t; } bool mitk::VirtualTrackingDevice::StartTracking() { if (this->GetState() != Ready) return false; this->SetState(Tracking); // go to mode Tracking this->m_StopTrackingMutex.lock(); this->m_StopTracking = false; this->m_StopTrackingMutex.unlock(); mitk::IGTTimeStamp::GetInstance()->Start(this); if (m_Thread.joinable()) m_Thread.detach(); m_Thread = std::thread(&VirtualTrackingDevice::ThreadStartTracking, this); // start a new thread that executes the TrackTools() method return true; } bool mitk::VirtualTrackingDevice::StopTracking() { if (this->GetState() == Tracking) // Only if the object is in the correct state { m_StopTrackingMutex.lock(); // m_StopTracking is used by two threads, so we have to ensure correct thread handling m_StopTracking = true; m_StopTrackingMutex.unlock(); m_TrackingFinishedMutex.lock(); this->SetState(Ready); m_TrackingFinishedMutex.unlock(); } mitk::IGTTimeStamp::GetInstance()->Stop(this); return true; } unsigned int mitk::VirtualTrackingDevice::GetToolCount() const { std::lock_guard lock(m_ToolsMutex); // lock and unlock the mutex return static_cast(this->m_AllTools.size()); } mitk::TrackingTool* mitk::VirtualTrackingDevice::GetTool(unsigned int toolNumber) const { std::lock_guard lock(m_ToolsMutex); // lock and unlock the mutex if (toolNumber < m_AllTools.size()) return this->m_AllTools.at(toolNumber); return nullptr; } bool mitk::VirtualTrackingDevice::OpenConnection() { if (m_NumberOfControlPoints < 1) { mitkThrowException(mitk::IGTException) << "to few control points for spline interpolation"; } srand(time(nullptr)); //Init random number generator this->SetState(Ready); return true; } void mitk::VirtualTrackingDevice::InitializeSpline(mitk::VirtualTrackingTool* t) { if (t == nullptr) return; typedef mitk::VirtualTrackingTool::SplineType SplineType; /* create random control points */ SplineType::ControlPointListType controlPoints; controlPoints.reserve(m_NumberOfControlPoints + 1); controlPoints.push_back(this->GetRandomPoint()); // insert point 0 double length = 0.0; // estimate spline length by calculating line segments lengths for (unsigned int i = 1; i < m_NumberOfControlPoints - 1; ++i) // set points 1..n-2 { SplineType::ControlPointType pos; pos = this->GetRandomPoint(); length += controlPoints.at(i - 1).EuclideanDistanceTo(pos); controlPoints.push_back(pos); } controlPoints.push_back(controlPoints.at(0)); // close spline --> insert point last control point with same value as first control point length += controlPoints.at(controlPoints.size() - 2).EuclideanDistanceTo(controlPoints.at(controlPoints.size() - 1)); /* Create knot list. TODO: rethink knot list values and list size. Is there a better solution? */ SplineType::KnotListType knotList; knotList.push_back(0.0); for (unsigned int i = 1; i < controlPoints.size() + t->GetSpline()->GetSplineOrder() + 1; ++i) knotList.push_back(i); knotList.push_back(controlPoints.size() + t->GetSpline()->GetSplineOrder() + 1); t->GetSpline()->SetControlPoints(controlPoints); t->GetSpline()->SetKnots(knotList); t->SetSplineLength(length); } bool mitk::VirtualTrackingDevice::CloseConnection() { bool returnValue = true; if (this->GetState() == Setup) return true; this->SetState(Setup); return returnValue; } mitk::ScalarType mitk::VirtualTrackingDevice::GetSplineChordLength(unsigned int idx) { mitk::VirtualTrackingTool* t = this->GetInternalTool(idx); if (t != nullptr) return t->GetSplineLength(); else throw std::invalid_argument("invalid index"); } void mitk::VirtualTrackingDevice::SetToolSpeed(unsigned int idx, mitk::ScalarType roundsPerSecond) { if (roundsPerSecond < 0.0001) throw std::invalid_argument("Minimum tool speed is 0.0001 rounds per second"); mitk::VirtualTrackingTool* t = this->GetInternalTool(idx); if (t != nullptr) t->SetVelocity(roundsPerSecond); else throw std::invalid_argument("invalid index"); } mitk::VirtualTrackingTool* mitk::VirtualTrackingDevice::GetInternalTool(unsigned int idx) { std::lock_guard toolsMutexLockHolder(m_ToolsMutex); // lock and unlock the mutex if (idx < m_AllTools.size()) return m_AllTools.at(idx); else return nullptr; } void mitk::VirtualTrackingDevice::TrackTools() { /* lock the TrackingFinishedMutex to signal that the execution rights are now transfered to the tracking thread */ std::lock_guard trackingFinishedLockHolder(m_TrackingFinishedMutex); // keep lock until end of scope if (this->GetState() != Tracking) return; bool localStopTracking; // Because m_StopTracking is used by two threads, access has to be guarded by a mutex. To minimize thread locking, a local copy is used here this->m_StopTrackingMutex.lock(); // update the local copy of m_StopTracking localStopTracking = this->m_StopTracking; this->m_StopTrackingMutex.unlock(); mitk::ScalarType t = 0.0; while ((this->GetState() == Tracking) && (localStopTracking == false)) { //for (ToolContainer::iterator itAllTools = m_AllTools.begin(); itAllTools != m_AllTools.end(); itAllTools++) for (unsigned int i = 0; i < this->GetToolCount(); ++i) // use mutexed methods to access tool container { mitk::VirtualTrackingTool::Pointer currentTool = this->GetInternalTool(i); mitk::VirtualTrackingTool::SplineType::PointType pos; /* calculate tool position with spline interpolation */ pos = currentTool->GetSpline()->EvaluateSpline(t); mitk::Point3D mp; mitk::itk2vtk(pos, mp); // convert from SplineType::PointType to mitk::Point3D //Add Gaussian Noise to Tracking Coordinates if enabled if (this->m_GaussianNoiseEnabled) { std::random_device rd; std::mt19937 generator(rd()); std::normal_distribution dist(this->m_MeanDistributionParam, this->m_DeviationDistributionParam); double noise = dist(generator); - mp = mp + noise; + mp = mp + mitk::Vector(noise); } currentTool->SetPosition(mp); // Currently, a constant speed is used. TODO: use tool velocity setting t += 0.001; if (t >= 1.0) t = 0.0; mitk::Quaternion quat; /* fix quaternion rotation */ quat.x() = 0.0; quat.y() = 0.0; quat.z() = 0.0; quat.r() = 1.0; quat.normalize(); currentTool->SetOrientation(quat); // TODO: rotate once per cycle around a fixed rotation vector currentTool->SetTrackingError(2 * (rand() / (RAND_MAX + 1.0))); // tracking error in 0 .. 2 Range currentTool->SetDataValid(true); currentTool->Modified(); } itksys::SystemTools::Delay(m_RefreshRate); /* Update the local copy of m_StopTracking */ this->m_StopTrackingMutex.lock(); localStopTracking = m_StopTracking; this->m_StopTrackingMutex.unlock(); } // tracking ends if we pass this line } void mitk::VirtualTrackingDevice::ThreadStartTracking() { this->TrackTools(); } mitk::VirtualTrackingDevice::ControlPointType mitk::VirtualTrackingDevice::GetRandomPoint() { ControlPointType pos; pos[0] = m_Bounds[0] + (m_Bounds[1] - m_Bounds[0]) * (rand() / (RAND_MAX + 1.0)); // X = xMin + xRange * (random number between 0 and 1) pos[1] = m_Bounds[2] + (m_Bounds[3] - m_Bounds[2]) * (rand() / (RAND_MAX + 1.0)); // Y pos[2] = m_Bounds[4] + (m_Bounds[5] - m_Bounds[4]) * (rand() / (RAND_MAX + 1.0)); // Z return pos; } void mitk::VirtualTrackingDevice::EnableGaussianNoise() { this->m_GaussianNoiseEnabled = true; } void mitk::VirtualTrackingDevice::DisableGaussianNoise() { this->m_GaussianNoiseEnabled = false; } void mitk::VirtualTrackingDevice::SetParamsForGaussianNoise(double meanDistribution, double deviationDistribution) { this->m_MeanDistributionParam = meanDistribution; this->m_DeviationDistributionParam = deviationDistribution; } double mitk::VirtualTrackingDevice::GetDeviationDistribution() { return m_DeviationDistributionParam; } double mitk::VirtualTrackingDevice::GetMeanDistribution() { return m_MeanDistributionParam; } diff --git a/Modules/ImageStatistics/CMakeLists.txt b/Modules/ImageStatistics/CMakeLists.txt index 1cf692bc6f..8647f7a8f0 100644 --- a/Modules/ImageStatistics/CMakeLists.txt +++ b/Modules/ImageStatistics/CMakeLists.txt @@ -1,8 +1,8 @@ mitk_create_module( DEPENDS MitkImageExtraction MitkPlanarFigure MitkMultilabel - PACKAGE_DEPENDS PRIVATE ITK|VTK VTK|IOImage + PACKAGE_DEPENDS PRIVATE ITK|Convolution+VTK VTK|IOImage ) if(BUILD_TESTING) add_subdirectory(Testing) endif() diff --git a/Modules/ImageStatistics/mitkMaskUtilities.tpp b/Modules/ImageStatistics/mitkMaskUtilities.tpp index f1d52ea2aa..ab740cb9f7 100644 --- a/Modules/ImageStatistics/mitkMaskUtilities.tpp +++ b/Modules/ImageStatistics/mitkMaskUtilities.tpp @@ -1,194 +1,194 @@ /*============================================================================ 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 MITKMASKUTIL_TPP #define MITKMASKUTIL_TPP #include #include #include #include #include namespace mitk { template void MaskUtilities::SetImage(const ImageType* image) { if (image != m_Image) { m_Image = image; } } template void MaskUtilities::SetMask(const MaskType* mask) { if (mask != m_Mask) { m_Mask = mask; } } template bool MaskUtilities::CheckMaskSanity() { if (m_Mask==nullptr || m_Image==nullptr) { MITK_ERROR << "Set an image and a mask first"; } typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::PointType PointType; typedef typename ImageType::DirectionType DirectionType; bool maskSanity = true; if (m_Mask==nullptr) { MITK_ERROR << "Something went wrong when casting the mitk mask image to an itk mask image. Do the mask and the input image have the same dimension?"; // note to self: We could try to convert say a 2d mask to a 3d mask if the image is 3d. (mask and image dimension have to match.) } // check direction DirectionType imageDirection = m_Image->GetDirection(); DirectionType maskDirection = m_Mask->GetDirection(); for(unsigned int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for(unsigned int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if (fabs(differenceDirection) > MASK_SUITABILITY_TOLERANCE_DIRECTION) { maskSanity = false; MITK_INFO << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")"; } } } // check spacing - PointType imageSpacing = m_Image->GetSpacing(); - PointType maskSpacing = m_Mask->GetSpacing(); + PointType imageSpacing(m_Image->GetSpacing().GetDataPointer()); + PointType maskSpacing(m_Mask->GetSpacing().GetDataPointer()); for (unsigned int i = 0; i < VImageDimension; i++) { if ( fabs( maskSpacing[i] - imageSpacing[i] ) > MASK_SUITABILITY_TOLERANCE_COORDINATE ) { maskSanity = false; MITK_INFO << "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing; } } // check alignment // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = m_Image->GetOrigin(); PointType maskOrigin = m_Mask->GetOrigin(); typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; m_Image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); m_Image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); // misalignment must be a multiple (int) of spacing in that direction if ( fmod(misalignment,imageSpacing[i]) > MASK_SUITABILITY_TOLERANCE_COORDINATE) { maskSanity = false; MITK_INFO << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << fmod(misalignment,imageSpacing[i]) << ")"; } } // mask must be completely inside image region // Make sure that mask region is contained within image region if ( m_Mask!=nullptr && !m_Image->GetLargestPossibleRegion().IsInside( m_Mask->GetLargestPossibleRegion() ) ) { maskSanity = false; MITK_INFO << "Mask region needs to be inside of image region! (Image region: " << m_Image->GetLargestPossibleRegion() << "; Mask region: " << m_Mask->GetLargestPossibleRegion() << ")"; } return maskSanity; } template typename MaskUtilities::ImageType::ConstPointer MaskUtilities::ExtractMaskImageRegion() { if (m_Mask==nullptr || m_Image==nullptr) { MITK_ERROR << "Set an image and a mask first"; } bool maskSanity = CheckMaskSanity(); if (!maskSanity) { MITK_ERROR << "Mask and image are not compatible"; } typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; typename ImageType::SizeType imageSize = m_Image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = m_Mask->GetBufferedRegion().GetSize(); typename itk::Image::ConstPointer resultImg; bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); typename MaskType::PointType maskOrigin = m_Mask->GetOrigin(); typename ImageType::PointType imageOrigin = m_Image->GetOrigin(); typename MaskType::SpacingType maskSpacing = m_Mask->GetSpacing(); typename ImageType::RegionType extractionRegion; typename ImageType::IndexType extractionRegionIndex; for (unsigned int i=0; i < maskOrigin.GetPointDimension(); i++) { extractionRegionIndex[i] = (maskOrigin[i] - imageOrigin[i]) / maskSpacing[i]; } extractionRegion.SetIndex(extractionRegionIndex); extractionRegion.SetSize(m_Mask->GetLargestPossibleRegion().GetSize()); extractImageFilter->SetInput( m_Image ); extractImageFilter->SetExtractionRegion( extractionRegion ); extractImageFilter->SetCoordinateTolerance(MASK_SUITABILITY_TOLERANCE_COORDINATE); extractImageFilter->SetDirectionTolerance(MASK_SUITABILITY_TOLERANCE_DIRECTION); extractImageFilter->Update(); auto extractedImg = extractImageFilter->GetOutput(); extractedImg->SetOrigin(m_Mask->GetOrigin()); extractedImg->SetLargestPossibleRegion(m_Mask->GetLargestPossibleRegion()); extractedImg->SetBufferedRegion(m_Mask->GetBufferedRegion()); resultImg = extractedImg; } else { resultImg = m_Image; } return resultImg; } } #endif diff --git a/Modules/Segmentation/Testing/mitkOverwriteSliceFilterObliquePlaneTest.cpp b/Modules/Segmentation/Testing/mitkOverwriteSliceFilterObliquePlaneTest.cpp index a6bb783040..e5839fa022 100644 --- a/Modules/Segmentation/Testing/mitkOverwriteSliceFilterObliquePlaneTest.cpp +++ b/Modules/Segmentation/Testing/mitkOverwriteSliceFilterObliquePlaneTest.cpp @@ -1,254 +1,254 @@ /*============================================================================ 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 #include #include #include #include #include #include #include #include #include #include #include int ObliquePlaneTestVolumeSize = 128; static void OverwriteObliquePlaneTest(mitk::Image *workingImage, mitk::Image *refImg) { /*==============TEST WITHOUT MITK CONVERTION=============================*/ /* ============= setup plane ============*/ auto sliceindex = (int)(ObliquePlaneTestVolumeSize / 2); // rand() % 32; bool isFrontside = true; bool isRotated = false; mitk::PlaneGeometry::Pointer obliquePlane = mitk::PlaneGeometry::New(); obliquePlane->InitializeStandardPlane( workingImage->GetGeometry(), mitk::AnatomicalPlane::Axial, sliceindex, isFrontside, isRotated); mitk::Point3D origin = obliquePlane->GetOrigin(); mitk::Vector3D normal; normal = obliquePlane->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 obliquePlane->SetOrigin(origin); mitk::Vector3D rotationVector = obliquePlane->GetAxisVector(0); rotationVector.Normalize(); float degree = 45.0; auto op = new mitk::RotationOperation(mitk::OpROTATE, obliquePlane->GetCenter(), rotationVector, degree); obliquePlane->ExecuteOperation(op); delete op; /* ============= extract slice ============*/ mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(); slicer->SetInput(workingImage); slicer->SetWorldGeometry(obliquePlane); slicer->SetVtkOutputRequest(true); slicer->Modified(); slicer->Update(); vtkSmartPointer slice = vtkSmartPointer::New(); slice = slicer->GetVtkOutput(); /* ============= overwrite slice ============*/ vtkSmartPointer resliceIdx = vtkSmartPointer::New(); mitk::ExtractSliceFilter::Pointer overwriter = mitk::ExtractSliceFilter::New(resliceIdx); resliceIdx->SetOverwriteMode(true); resliceIdx->SetInputSlice(slice); resliceIdx->Modified(); overwriter->SetInput(workingImage); overwriter->SetWorldGeometry(obliquePlane); overwriter->SetVtkOutputRequest(true); overwriter->Modified(); overwriter->Update(); typedef mitk::ImagePixelReadAccessor ReadAccessorType; ReadAccessorType refImgReadAccessor(refImg); ReadAccessorType workingImgReadAccessor(workingImage); /* ============= check ref == working ============*/ bool areSame = true; itk::Index<3> id; id[0] = id[1] = id[2] = 0; for (int x = 0; x < ObliquePlaneTestVolumeSize; ++x) { id[0] = x; for (int y = 0; y < ObliquePlaneTestVolumeSize; ++y) { id[1] = y; for (int z = 0; z < ObliquePlaneTestVolumeSize; ++z) { id[2] = z; areSame = refImgReadAccessor.GetPixelByIndex(id) == workingImgReadAccessor.GetPixelByIndex(id); if (!areSame) goto stop; } } } stop: MITK_TEST_CONDITION(areSame, "comparing images (no mitk convertion) [oblique]"); /*==============TEST WITH MITK CONVERTION=============================*/ /* ============= extract slice ============*/ mitk::ExtractSliceFilter::Pointer slicer2 = mitk::ExtractSliceFilter::New(); slicer2->SetInput(workingImage); slicer2->SetWorldGeometry(obliquePlane); slicer2->Modified(); slicer2->Update(); mitk::Image::Pointer sliceInMitk = slicer2->GetOutput(); vtkSmartPointer slice2 = vtkSmartPointer::New(); slice2 = sliceInMitk->GetVtkImageData(); /* ============= overwrite slice ============*/ vtkSmartPointer resliceIdx2 = vtkSmartPointer::New(); mitk::ExtractSliceFilter::Pointer overwriter2 = mitk::ExtractSliceFilter::New(resliceIdx2); resliceIdx2->SetOverwriteMode(true); resliceIdx2->SetInputSlice(slice2); resliceIdx2->Modified(); overwriter2->SetInput(workingImage); overwriter2->SetWorldGeometry(obliquePlane); overwriter2->SetVtkOutputRequest(true); overwriter2->Modified(); overwriter2->Update(); /* ============= check ref == working ============*/ areSame = true; id[0] = id[1] = id[2] = 0; for (int x = 0; x < ObliquePlaneTestVolumeSize; ++x) { id[0] = x; for (int y = 0; y < ObliquePlaneTestVolumeSize; ++y) { id[1] = y; for (int z = 0; z < ObliquePlaneTestVolumeSize; ++z) { id[2] = z; areSame = refImgReadAccessor.GetPixelByIndex(id) == workingImgReadAccessor.GetPixelByIndex(id); if (!areSame) goto stop2; } } } stop2: MITK_TEST_CONDITION(areSame, "comparing images (with mitk convertion) [oblique]"); /*==============TEST EDIT WITHOUT MITK CONVERTION=============================*/ /* ============= edit slice ============*/ int idX = std::abs(ObliquePlaneTestVolumeSize - 59); int idY = std::abs(ObliquePlaneTestVolumeSize - 23); int idZ = 0; int component = 0; double val = 33.0; slice->SetScalarComponentFromDouble(idX, idY, idZ, component, val); mitk::Vector3D indx; indx[0] = idX; indx[1] = idY; indx[2] = idZ; sliceInMitk->GetGeometry()->IndexToWorld(indx, indx); /* ============= overwrite slice ============*/ vtkSmartPointer resliceIdx3 = vtkSmartPointer::New(); resliceIdx3->SetOverwriteMode(true); resliceIdx3->SetInputSlice(slice); mitk::ExtractSliceFilter::Pointer overwriter3 = mitk::ExtractSliceFilter::New(resliceIdx3); overwriter3->SetInput(workingImage); overwriter3->SetWorldGeometry(obliquePlane); overwriter3->SetVtkOutputRequest(true); overwriter3->Modified(); overwriter3->Update(); /* ============= check ============*/ areSame = true; int x = 0; int y = 0; int z = 0; for (x = 0; x < ObliquePlaneTestVolumeSize; ++x) { id[0] = x; for (y = 0; y < ObliquePlaneTestVolumeSize; ++y) { id[1] = y; for (z = 0; z < ObliquePlaneTestVolumeSize; ++z) { id[2] = z; areSame = refImgReadAccessor.GetPixelByIndex(id) == workingImgReadAccessor.GetPixelByIndex(id); if (!areSame) goto stop3; } } } stop3: MITK_TEST_CONDITION(x == idX && y == z, "overwrited the right index [oblique]"); } /*================ #BEGIN test main ================*/ int mitkOverwriteSliceFilterObliquePlaneTest(int, char *[]) { MITK_TEST_BEGIN("mitkOverwriteSliceFilterObliquePlaneTest") typedef itk::Image ImageType; typedef itk::ImageRegionConstIterator ImageIterator; ImageType::Pointer image = ImageType::New(); ImageType::IndexType start; start[0] = start[1] = start[2] = 0; ImageType::SizeType size; size[0] = size[1] = size[2] = ObliquePlaneTestVolumeSize; ImageType::RegionType imgRegion; imgRegion.SetSize(size); imgRegion.SetIndex(start); image->SetRegions(imgRegion); - image->SetSpacing(1.0); + image->SetSpacing(mitk::Vector(1.0)); image->Allocate(); ImageIterator imageIterator(image, image->GetLargestPossibleRegion()); imageIterator.GoToBegin(); unsigned short pixelValue = 0; // fill the image with distinct values while (!imageIterator.IsAtEnd()) { image->SetPixel(imageIterator.GetIndex(), pixelValue); ++imageIterator; ++pixelValue; } /* end setup itk image */ mitk::Image::Pointer refImage; CastToMitkImage(image, refImage); mitk::Image::Pointer workingImg; CastToMitkImage(image, workingImg); OverwriteObliquePlaneTest(workingImg, refImage); MITK_TEST_END() } diff --git a/Modules/Segmentation/Testing/mitkOverwriteSliceFilterTest.cpp b/Modules/Segmentation/Testing/mitkOverwriteSliceFilterTest.cpp index 6b67791740..f522844be6 100644 --- a/Modules/Segmentation/Testing/mitkOverwriteSliceFilterTest.cpp +++ b/Modules/Segmentation/Testing/mitkOverwriteSliceFilterTest.cpp @@ -1,181 +1,181 @@ /*============================================================================ 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 #include #include #include #include #include #include #include #include #include #include #include int VolumeSize = 128; /*================ #BEGIN test main ================*/ int mitkOverwriteSliceFilterTest(int, char *[]) { MITK_TEST_BEGIN("mitkOverwriteSliceFilterTest") typedef itk::Image ImageType; typedef itk::ImageRegionConstIterator ImageIterator; ImageType::Pointer image = ImageType::New(); ImageType::IndexType start; start[0] = start[1] = start[2] = 0; ImageType::SizeType size; size[0] = size[1] = size[2] = VolumeSize; ImageType::RegionType imgRegion; imgRegion.SetSize(size); imgRegion.SetIndex(start); image->SetRegions(imgRegion); - image->SetSpacing(1.0); + image->SetSpacing(mitk::Vector(1.0)); image->Allocate(); ImageIterator imageIterator(image, image->GetLargestPossibleRegion()); imageIterator.GoToBegin(); unsigned short pixelValue = 0; // fill the image with distinct values while (!imageIterator.IsAtEnd()) { image->SetPixel(imageIterator.GetIndex(), pixelValue); ++imageIterator; ++pixelValue; } /* end setup itk image */ mitk::Image::Pointer referenceImage; CastToMitkImage(image, referenceImage); mitk::Image::Pointer workingImage; CastToMitkImage(image, workingImage); typedef mitk::ImagePixelReadAccessor ReadAccessorType; ReadAccessorType refImgReadAccessor(referenceImage); ReadAccessorType workingImgReadAccessor(workingImage); /* ============= setup plane ============*/ int sliceindex = 55; // rand() % 32; bool isFrontside = true; bool isRotated = false; mitk::PlaneGeometry::Pointer plane = mitk::PlaneGeometry::New(); plane->InitializeStandardPlane( workingImage->GetGeometry(), mitk::AnatomicalPlane::Axial, sliceindex, isFrontside, isRotated); mitk::Point3D origin = plane->GetOrigin(); mitk::Vector3D normal; normal = plane->GetNormal(); normal.Normalize(); origin += normal * 0.5; // pixelspacing is 1, so half the spacing is 0.5 plane->SetOrigin(origin); /* ============= extract slice ============*/ vtkSmartPointer resliceIdx = vtkSmartPointer::New(); mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(resliceIdx); slicer->SetInput(workingImage); slicer->SetWorldGeometry(plane); slicer->SetVtkOutputRequest(true); slicer->Modified(); slicer->Update(); vtkSmartPointer slice = vtkSmartPointer::New(); slice = slicer->GetVtkOutput(); /* ============= overwrite slice ============*/ resliceIdx->SetOverwriteMode(true); resliceIdx->Modified(); slicer->Modified(); slicer->Update(); // implicit overwrite /* ============= check ref == working ============*/ bool areSame = true; itk::Index<3> id; id[0] = id[1] = id[2] = 0; for (int x = 0; x < VolumeSize; ++x) { id[0] = x; for (int y = 0; y < VolumeSize; ++y) { id[1] = y; for (int z = 0; z < VolumeSize; ++z) { id[2] = z; areSame = refImgReadAccessor.GetPixelByIndex(id) == workingImgReadAccessor.GetPixelByIndex(id); if (!areSame) goto stop; } } } stop: MITK_TEST_CONDITION(areSame, "test overwrite unmodified slice"); /* ============= edit slice ============*/ int idX = std::abs(VolumeSize - 59); int idY = std::abs(VolumeSize - 23); int idZ = 0; int component = 0; double val = 33.0; slice->SetScalarComponentFromDouble(idX, idY, idZ, component, val); /* ============= overwrite slice ============*/ vtkSmartPointer resliceIdx2 = vtkSmartPointer::New(); resliceIdx2->SetOverwriteMode(true); resliceIdx2->SetInputSlice(slice); mitk::ExtractSliceFilter::Pointer slicer2 = mitk::ExtractSliceFilter::New(resliceIdx2); slicer2->SetInput(workingImage); slicer2->SetWorldGeometry(plane); slicer2->SetVtkOutputRequest(true); slicer2->Modified(); slicer2->Update(); /* ============= check ============*/ areSame = true; int xx = 0; int yy = 0; int zz = 0; for (xx = 0; xx < VolumeSize; ++xx) { id[0] = xx; for (yy = 0; yy < VolumeSize; ++yy) { id[1] = yy; for (zz = 0; zz < VolumeSize; ++zz) { id[2] = zz; areSame = refImgReadAccessor.GetPixelByIndex(id) == workingImgReadAccessor.GetPixelByIndex(id); if (!areSame) goto stop2; } } } stop2: // MITK_INFO << "index: [" << x << ", " << y << ", " << z << "]"; MITK_TEST_CONDITION(xx == idX && yy == idY && zz == sliceindex, "test overwrite modified slice"); MITK_TEST_END() } diff --git a/Modules/SurfaceInterpolation/Testing/mitkSurfaceInterpolationControllerTest.cpp b/Modules/SurfaceInterpolation/Testing/mitkSurfaceInterpolationControllerTest.cpp index 20faa7a6e4..26418f9785 100644 --- a/Modules/SurfaceInterpolation/Testing/mitkSurfaceInterpolationControllerTest.cpp +++ b/Modules/SurfaceInterpolation/Testing/mitkSurfaceInterpolationControllerTest.cpp @@ -1,1097 +1,1097 @@ /*============================================================================ 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 #include #include #include #include #include #include #include #include #include class mitkSurfaceInterpolationControllerTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkSurfaceInterpolationControllerTestSuite); MITK_TEST(TestSingleton); MITK_TEST(TestSetCurrentInterpolationSession); MITK_TEST(TestReplaceInterpolationSession); MITK_TEST(TestRemoveAllInterpolationSessions); MITK_TEST(TestRemoveInterpolationSession); MITK_TEST(TestOnSegmentationDeleted); MITK_TEST(TestSetCurrentInterpolationSession4D); MITK_TEST(TestReplaceInterpolationSession4D); MITK_TEST(TestRemoveAllInterpolationSessions4D); MITK_TEST(TestRemoveInterpolationSession4D); MITK_TEST(TestOnSegmentationDeleted4D); /// \todo Workaround for memory leak in TestAddNewContour. Bug 18096. vtkDebugLeaks::SetExitError(0); MITK_TEST(TestAddNewContour); MITK_TEST(TestRemoveContour); CPPUNIT_TEST_SUITE_END(); private: mitk::SurfaceInterpolationController::Pointer m_Controller; public: mitk::Image::Pointer createImage(unsigned int *dimensions) { mitk::Image::Pointer newImage = mitk::Image::New(); // mitk::LabelSetImage::Pointer newImage = mitk::LabelSetImage::New(); mitk::PixelType p_type = mitk::MakeScalarPixelType(); newImage->Initialize(p_type, 3, dimensions); return newImage; } mitk::LabelSetImage::Pointer createLabelSetImage(unsigned int *dimensions) { mitk::Image::Pointer image = createImage(dimensions); mitk::LabelSetImage::Pointer newImage = mitk::LabelSetImage::New(); newImage->InitializeByLabeledImage(image); return newImage; } mitk::Image::Pointer createImage4D(unsigned int *dimensions) { mitk::Image::Pointer newImage = mitk::Image::New(); mitk::PixelType p_type = mitk::MakeScalarPixelType(); newImage->Initialize(p_type, 4, dimensions); return newImage; } mitk::LabelSetImage::Pointer createLabelSetImage4D(unsigned int *dimensions) { mitk::Image::Pointer image = createImage4D(dimensions); mitk::LabelSetImage::Pointer newImage = mitk::LabelSetImage::New(); newImage->InitializeByLabeledImage(image); return newImage; } void setUp() override { m_Controller = mitk::SurfaceInterpolationController::GetInstance(); m_Controller->RemoveAllInterpolationSessions(); m_Controller->SetCurrentTimePoint(0); vtkSmartPointer polygonSource = vtkSmartPointer::New(); polygonSource->SetRadius(100); polygonSource->SetNumberOfSides(7); polygonSource->Update(); mitk::Surface::Pointer surface = mitk::Surface::New(); surface->SetVtkPolyData(polygonSource->GetOutput()); } void TestSingleton() { mitk::SurfaceInterpolationController::Pointer controller2 = mitk::SurfaceInterpolationController::GetInstance(); CPPUNIT_ASSERT_MESSAGE("SurfaceInterpolationController pointers are not equal!", m_Controller.GetPointer() == controller2.GetPointer()); } void TestSetCurrentInterpolationSession() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage(dimensions1); unsigned int dimensions2[] = {20, 10, 30}; mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage(dimensions2); // Test 1 m_Controller->SetCurrentInterpolationSession(segmentation_1); MITK_ASSERT_EQUAL( m_Controller->GetCurrentSegmentation(), segmentation_1->Clone(), "Segmentation images are not equal"); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_1.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); // Test 2 m_Controller->SetCurrentInterpolationSession(segmentation_2); MITK_ASSERT_EQUAL( m_Controller->GetCurrentSegmentation(), segmentation_2->Clone(), "Segmentation images are not equal"); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_2.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test 3 m_Controller->SetCurrentInterpolationSession(segmentation_1); MITK_ASSERT_EQUAL( m_Controller->GetCurrentSegmentation(), segmentation_1->Clone(), "Segmentation images are not equal"); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_1.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test 4 m_Controller->SetCurrentInterpolationSession(segmentation_1); MITK_ASSERT_EQUAL( m_Controller->GetCurrentSegmentation(), segmentation_1->Clone(), "Segmentation images are not equal"); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_1.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // // Test 5 m_Controller->SetCurrentInterpolationSession(nullptr); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().IsNull()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); } mitk::PlaneGeometry::Pointer GetPlaneGeometry() { mitk::Point3D origin; mitk::Vector3D right, bottom, normal, spacing; mitk::ScalarType width, height; mitk::ScalarType widthInMM, heightInMM, thicknessInMM; auto planegeometry = mitk::PlaneGeometry::New(); width = 100; widthInMM = width; height = 200; heightInMM = height; thicknessInMM = 1.0; mitk::FillVector3D(origin, 4.5, 7.3, 11.2); mitk::FillVector3D(right, widthInMM, 0, 0); mitk::FillVector3D(bottom, 0, heightInMM, 0); mitk::FillVector3D(normal, 0, 0, thicknessInMM); mitk::FillVector3D(spacing, 1.0, 1.0, thicknessInMM); planegeometry->InitializeStandardPlane(right, bottom); planegeometry->SetOrigin(origin); planegeometry->SetSpacing(spacing); return planegeometry; } void TestReplaceInterpolationSession() { // Create segmentation image unsigned int dimensions1[] = {10, 10, 10}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage(dimensions1); m_Controller->SetCurrentInterpolationSession(segmentation_1); // Create some contours double center_1[3] = {1.25f, 3.43f, 4.44f}; double normal_1[3] = {0.25f, 1.76f, 0.93f}; vtkSmartPointer p_source = vtkSmartPointer::New(); p_source->SetNumberOfSides(20); p_source->SetCenter(center_1); p_source->SetRadius(4); p_source->SetNormal(normal_1); p_source->Update(); vtkPolyData *poly_1 = p_source->GetOutput(); mitk::Surface::Pointer surf_1 = mitk::Surface::New(); surf_1->SetVtkPolyData(poly_1); vtkSmartPointer int1Array = vtkSmartPointer::New(); int1Array->InsertNextValue(1); int1Array->InsertNextValue(0); int1Array->InsertNextValue(0); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(int1Array); vtkSmartPointer double1Array = vtkSmartPointer::New(); double1Array->InsertNextValue(center_1[0]); double1Array->InsertNextValue(center_1[1]); double1Array->InsertNextValue(center_1[2]); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(double1Array); double center_2[3] = {4.0f, 4.0f, 4.0f}; double normal_2[3] = {1.0f, 0.0f, 0.0f}; vtkSmartPointer p_source_2 = vtkSmartPointer::New(); p_source_2->SetNumberOfSides(80); p_source_2->SetCenter(center_2); p_source_2->SetRadius(4); p_source_2->SetNormal(normal_2); p_source_2->Update(); vtkPolyData *poly_2 = p_source_2->GetOutput(); mitk::Surface::Pointer surf_2 = mitk::Surface::New(); surf_2->SetVtkPolyData(poly_2); vtkSmartPointer int2Array = vtkSmartPointer::New(); int2Array->InsertNextValue(1); int2Array->InsertNextValue(0); int2Array->InsertNextValue(0); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(int2Array); vtkSmartPointer doubleArray = vtkSmartPointer::New(); doubleArray->InsertNextValue(center_2[0]); doubleArray->InsertNextValue(center_2[1]); doubleArray->InsertNextValue(center_2[2]); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(doubleArray); std::vector surfaces; surfaces.push_back(surf_1); surfaces.push_back(surf_2); const mitk::PlaneGeometry * planeGeometry1 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry2 = GetPlaneGeometry(); std::vector planeGeometries; planeGeometries.push_back(planeGeometry1); planeGeometries.push_back(planeGeometry2); // Add contours m_Controller->AddNewContours(surfaces, planeGeometries, true); // Check if all contours are there mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo1; mitk::ScalarType n[3]; vtkPolygon::ComputeNormal(surf_1->GetVtkPolyData()->GetPoints(), n); contourInfo1.ContourNormal = n; contourInfo1.ContourPoint = center_1; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo2; vtkPolygon::ComputeNormal(surf_2->GetVtkPolyData()->GetPoints(), n); contourInfo2.ContourNormal = n; contourInfo2.ContourPoint = center_2; const mitk::Surface *contour_1 = m_Controller->GetContour(contourInfo1); const mitk::Surface *contour_2 = m_Controller->GetContour(contourInfo2); CPPUNIT_ASSERT_MESSAGE("Wrong number of contours!", m_Controller->GetNumberOfContours() == 2); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_1->GetVtkPolyData()), *(contour_1->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_2->GetVtkPolyData()), *(contour_2->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage(dimensions1); bool success = m_Controller->ReplaceInterpolationSession(segmentation_1, segmentation_2); CPPUNIT_ASSERT_MESSAGE("Replace session failed!", success == true); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_2.GetPointer()); unsigned int dimensions2[] = {10, 20, 10}; mitk::Image::Pointer segmentation_3 = createLabelSetImage(dimensions2); success = m_Controller->ReplaceInterpolationSession(segmentation_1, segmentation_3); CPPUNIT_ASSERT_MESSAGE("Replace session failed!", success == false); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_2.GetPointer()); } void TestRemoveAllInterpolationSessions() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10}; auto segmentation_1 = createLabelSetImage(dimensions1); unsigned int dimensions2[] = {20, 10, 30}; auto segmentation_2 = createLabelSetImage(dimensions2); // Test 1 m_Controller->SetCurrentInterpolationSession(segmentation_1); m_Controller->SetCurrentInterpolationSession(segmentation_2); m_Controller->RemoveAllInterpolationSessions(); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 0", m_Controller->GetNumberOfInterpolationSessions() == 0); } void TestRemoveInterpolationSession() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage(dimensions1); unsigned int dimensions2[] = {20, 10, 30}; mitk::Image::Pointer segmentation_2 = createLabelSetImage(dimensions2); // Test 1 m_Controller->SetCurrentInterpolationSession(segmentation_1); m_Controller->SetCurrentInterpolationSession(segmentation_2); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test current segmentation should not be null if another one was removed m_Controller->RemoveInterpolationSession(segmentation_1); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_2.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Current segmentation is null after another one was removed", m_Controller->GetCurrentSegmentation().IsNotNull()); m_Controller->SetCurrentInterpolationSession(segmentation_1); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test current segmentation should not be null if another one was removed m_Controller->RemoveInterpolationSession(segmentation_1); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); CPPUNIT_ASSERT_MESSAGE("Current segmentation is not null after session was removed", m_Controller->GetCurrentSegmentation().IsNull()); } void TestOnSegmentationDeleted() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10}; auto segmentation_1 = createLabelSetImage(dimensions1); m_Controller->SetCurrentInterpolationSession(segmentation_1); segmentation_1 = nullptr; CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 0", m_Controller->GetNumberOfInterpolationSessions() == 0); } void TestAddNewContour() { // Create segmentation image unsigned int dimensions1[] = {10, 10, 10}; mitk::Image::Pointer segmentation_1 = createLabelSetImage(dimensions1); m_Controller->SetCurrentInterpolationSession(segmentation_1); // Create some contours double center_1[3] = {1.25f, 3.43f, 4.44f}; double normal_1[3] = {0.25f, 1.76f, 0.93f}; vtkSmartPointer p_source = vtkSmartPointer::New(); p_source->SetNumberOfSides(20); p_source->SetCenter(center_1); p_source->SetRadius(4); p_source->SetNormal(normal_1); p_source->Update(); vtkPolyData *poly_1 = p_source->GetOutput(); mitk::Surface::Pointer surf_1 = mitk::Surface::New(); surf_1->SetVtkPolyData(poly_1); vtkSmartPointer int1Array = vtkSmartPointer::New(); int1Array->InsertNextValue(1); int1Array->InsertNextValue(0); int1Array->InsertNextValue(0); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(int1Array); vtkSmartPointer double1Array = vtkSmartPointer::New(); double1Array->InsertNextValue(center_1[0]); double1Array->InsertNextValue(center_1[1]); double1Array->InsertNextValue(center_1[2]); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(double1Array); double center_2[3] = {4.0f, 4.0f, 4.0f}; double normal_2[3] = {1.0f, 0.0f, 0.0f}; vtkSmartPointer p_source_2 = vtkSmartPointer::New(); p_source_2->SetNumberOfSides(80); p_source_2->SetCenter(center_2); p_source_2->SetRadius(4); p_source_2->SetNormal(normal_2); p_source_2->Update(); vtkPolyData *poly_2 = p_source_2->GetOutput(); mitk::Surface::Pointer surf_2 = mitk::Surface::New(); surf_2->SetVtkPolyData(poly_2); vtkSmartPointer int2Array = vtkSmartPointer::New(); int2Array->InsertNextValue(1); int2Array->InsertNextValue(0); int2Array->InsertNextValue(0); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(int2Array); vtkSmartPointer double2Array = vtkSmartPointer::New(); double2Array->InsertNextValue(center_2[0]); double2Array->InsertNextValue(center_2[1]); double2Array->InsertNextValue(center_2[2]); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(double2Array); double center_3[3] = {4.0f, 4.0f, 3.0f}; double normal_3[3] = {0.0f, 0.0f, 1.0f}; vtkSmartPointer p_source_3 = vtkSmartPointer::New(); p_source_3->SetNumberOfSides(10); p_source_3->SetCenter(center_3); p_source_3->SetRadius(4); p_source_3->SetNormal(normal_3); p_source_3->Update(); vtkPolyData *poly_3 = p_source_3->GetOutput(); mitk::Surface::Pointer surf_3 = mitk::Surface::New(); surf_3->SetVtkPolyData(poly_3); vtkSmartPointer int3Array = vtkSmartPointer::New(); int3Array->InsertNextValue(1); int3Array->InsertNextValue(0); int3Array->InsertNextValue(0); surf_3->GetVtkPolyData()->GetFieldData()->AddArray(int3Array); vtkSmartPointer double3Array = vtkSmartPointer::New(); double3Array->InsertNextValue(center_3[0]); double3Array->InsertNextValue(center_3[1]); double3Array->InsertNextValue(center_3[2]); surf_3->GetVtkPolyData()->GetFieldData()->AddArray(double3Array); std::vector surfaces; surfaces.push_back(surf_1); surfaces.push_back(surf_2); surfaces.push_back(surf_3); const mitk::PlaneGeometry * planeGeometry1 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry2 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry3 = GetPlaneGeometry(); std::vector planeGeometries; planeGeometries.push_back(planeGeometry1); planeGeometries.push_back(planeGeometry2); planeGeometries.push_back(planeGeometry3); // Add contours m_Controller->AddNewContours(surfaces, planeGeometries, true); mitk::ScalarType n[3]; // Check if all contours are there mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo1; vtkPolygon::ComputeNormal(surf_1->GetVtkPolyData()->GetPoints(), n); contourInfo1.ContourNormal = n; contourInfo1.ContourPoint = center_1; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo2; vtkPolygon::ComputeNormal(surf_2->GetVtkPolyData()->GetPoints(), n); contourInfo2.ContourNormal = n; contourInfo2.ContourPoint = center_2; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo3; vtkPolygon::ComputeNormal(surf_3->GetVtkPolyData()->GetPoints(), n); contourInfo3.ContourNormal = n; contourInfo3.ContourPoint = center_3; const mitk::Surface *contour_1 = m_Controller->GetContour(contourInfo1); const mitk::Surface *contour_2 = m_Controller->GetContour(contourInfo2); const mitk::Surface *contour_3 = m_Controller->GetContour(contourInfo3); CPPUNIT_ASSERT_MESSAGE("Wrong number of contours!", m_Controller->GetNumberOfContours() == 3); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_1->GetVtkPolyData()), *(contour_1->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_2->GetVtkPolyData()), *(contour_2->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_3->GetVtkPolyData()), *(contour_3->GetVtkPolyData()), 0.000001, true)); // Create another segmentation image unsigned int dimensions2[] = {20, 20, 20}; mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage(dimensions2); m_Controller->SetCurrentInterpolationSession(segmentation_2); // Create some contours double center_4[3] = {10.0f, 10.0f, 10.0f}; double normal_4[3] = {0.0f, 1.0f, 0.0f}; vtkSmartPointer p_source_4 = vtkSmartPointer::New(); p_source_4->SetNumberOfSides(8); p_source_4->SetCenter(center_4); p_source_4->SetRadius(5); p_source_4->SetNormal(normal_4); p_source_4->Update(); vtkPolyData *poly_4 = p_source_4->GetOutput(); mitk::Surface::Pointer surf_4 = mitk::Surface::New(); surf_4->SetVtkPolyData(poly_4); vtkSmartPointer int4Array = vtkSmartPointer::New(); int4Array->InsertNextValue(2); int4Array->InsertNextValue(0); int4Array->InsertNextValue(0); surf_4->GetVtkPolyData()->GetFieldData()->AddArray(int4Array); vtkSmartPointer double4Array = vtkSmartPointer::New(); double4Array->InsertNextValue(center_4[0]); double4Array->InsertNextValue(center_4[1]); double4Array->InsertNextValue(center_4[2]); surf_4->GetVtkPolyData()->GetFieldData()->AddArray(double4Array); double center_5[3] = {3.0f, 10.0f, 10.0f}; double normal_5[3] = {1.0f, 0.0f, 0.0f}; vtkSmartPointer p_source_5 = vtkSmartPointer::New(); p_source_5->SetNumberOfSides(16); p_source_5->SetCenter(center_5); p_source_5->SetRadius(8); p_source_5->SetNormal(normal_5); p_source_5->Update(); vtkPolyData *poly_5 = p_source_5->GetOutput(); mitk::Surface::Pointer surf_5 = mitk::Surface::New(); surf_5->SetVtkPolyData(poly_5); vtkSmartPointer int5Array = vtkSmartPointer::New(); int5Array->InsertNextValue(2); int5Array->InsertNextValue(0); int5Array->InsertNextValue(0); surf_5->GetVtkPolyData()->GetFieldData()->AddArray(int5Array); vtkSmartPointer double5Array = vtkSmartPointer::New(); double5Array->InsertNextValue(center_5[0]); double5Array->InsertNextValue(center_5[1]); double5Array->InsertNextValue(center_5[2]); surf_5->GetVtkPolyData()->GetFieldData()->AddArray(double5Array); double center_6[3] = {10.0f, 10.0f, 3.0f}; double normal_6[3] = {0.0f, 0.0f, 1.0f}; vtkSmartPointer p_source_6 = vtkSmartPointer::New(); p_source_6->SetNumberOfSides(100); p_source_6->SetCenter(center_6); p_source_6->SetRadius(5); p_source_6->SetNormal(normal_6); p_source_6->Update(); vtkPolyData *poly_6 = p_source_6->GetOutput(); mitk::Surface::Pointer surf_6 = mitk::Surface::New(); surf_6->SetVtkPolyData(poly_6); vtkSmartPointer int6Array = vtkSmartPointer::New(); int6Array->InsertNextValue(2); int6Array->InsertNextValue(0); int6Array->InsertNextValue(0); surf_6->GetVtkPolyData()->GetFieldData()->AddArray(int6Array); vtkSmartPointer double6Array = vtkSmartPointer::New(); double6Array->InsertNextValue(center_6[0]); double6Array->InsertNextValue(center_6[1]); double6Array->InsertNextValue(center_6[2]); surf_6->GetVtkPolyData()->GetFieldData()->AddArray(double6Array); mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo4; vtkPolygon::ComputeNormal(surf_4->GetVtkPolyData()->GetPoints(), n); contourInfo4.ContourNormal = n; contourInfo4.ContourPoint = center_4; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo5; vtkPolygon::ComputeNormal(surf_5->GetVtkPolyData()->GetPoints(), n); contourInfo5.ContourNormal = n; contourInfo5.ContourPoint = center_5; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo6; vtkPolygon::ComputeNormal(surf_6->GetVtkPolyData()->GetPoints(), n); contourInfo6.ContourNormal = n; contourInfo6.ContourPoint = center_6; const mitk::PlaneGeometry * planeGeometry4 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry5 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry6 = GetPlaneGeometry(); std::vector surfaces2; surfaces2.push_back(surf_4); surfaces2.push_back(surf_5); surfaces2.push_back(surf_6); std::vector planeGeometries2; planeGeometries2.push_back(planeGeometry4); planeGeometries2.push_back(planeGeometry5); planeGeometries2.push_back(planeGeometry6); m_Controller->AddNewContours(surfaces2, planeGeometries2, true); // Check if all contours are there auto contour_4 = m_Controller->GetContour(contourInfo4); auto contour_5 = m_Controller->GetContour(contourInfo5); auto contour_6 = m_Controller->GetContour(contourInfo6); CPPUNIT_ASSERT_MESSAGE("Wrong number of contours!", m_Controller->GetNumberOfContours() == 3); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_4->GetVtkPolyData()), *(contour_4->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_5->GetVtkPolyData()), *(contour_5->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_6->GetVtkPolyData()), *(contour_6->GetVtkPolyData()), 0.000001, true)); // Modify some contours vtkSmartPointer p_source_7 = vtkSmartPointer::New(); p_source_7->SetNumberOfSides(200); p_source_7->SetCenter(3.0, 10.0, 10.0); p_source_7->SetRadius(5); p_source_7->SetNormal(1, 0, 0); p_source_7->Update(); vtkPolyData *poly_7 = p_source_7->GetOutput(); mitk::Surface::Pointer surf_7 = mitk::Surface::New(); surf_7->SetVtkPolyData(poly_7); vtkSmartPointer int7Array = vtkSmartPointer::New(); int7Array->InsertNextValue(2); int7Array->InsertNextValue(0); int7Array->InsertNextValue(0); surf_7->GetVtkPolyData()->GetFieldData()->AddArray(int7Array); vtkSmartPointer double7Array = vtkSmartPointer::New(); double7Array->InsertNextValue(3.0); double7Array->InsertNextValue(10.0); double7Array->InsertNextValue(10.0); surf_7->GetVtkPolyData()->GetFieldData()->AddArray(double7Array); std::vector surfaces3; surfaces3.push_back(surf_7); const mitk::PlaneGeometry * planeGeometry7 = GetPlaneGeometry(); std::vector planeGeometries3; planeGeometries3.push_back(planeGeometry7); m_Controller->AddNewContours(surfaces3, planeGeometries3, true); mitk::ScalarType center_7[3]; center_7[0] = 3.0; center_7[1] = 10.0; center_7[2] = 10.0; vtkPolygon::ComputeNormal(surf_7->GetVtkPolyData()->GetPoints(), n); contourInfo5.ContourNormal = n; contourInfo5.ContourPoint = center_7; auto contour_7 = m_Controller->GetContour(contourInfo5); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_7->GetVtkPolyData()), *(contour_7->GetVtkPolyData()), 0.000001, true)); // Change session and test if all contours are available m_Controller->SetCurrentInterpolationSession(segmentation_1); auto contour_8 = m_Controller->GetContour(contourInfo1); auto contour_9 = m_Controller->GetContour(contourInfo2); auto contour_10 = m_Controller->GetContour(contourInfo3); CPPUNIT_ASSERT_MESSAGE("Wrong number of contours!", m_Controller->GetNumberOfContours() == 3); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_1->GetVtkPolyData()), *(contour_8->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_2->GetVtkPolyData()), *(contour_9->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_3->GetVtkPolyData()), *(contour_10->GetVtkPolyData()), 0.000001, true)); } void TestRemoveContour() { // Create segmentation image unsigned int dimensions1[] = {12, 12, 12}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage(dimensions1); m_Controller->SetCurrentInterpolationSession(segmentation_1); // Create some contours double center_1[3] = {4.0f, 4.0f, 4.0f}; double normal_1[3] = {0.0f, 1.0f, 0.0f}; vtkSmartPointer p_source = vtkSmartPointer::New(); p_source->SetNumberOfSides(20); p_source->SetCenter(center_1); p_source->SetRadius(4); p_source->SetNormal(normal_1); p_source->Update(); vtkPolyData *poly_1 = p_source->GetOutput(); mitk::Surface::Pointer surf_1 = mitk::Surface::New(); surf_1->SetVtkPolyData(poly_1); vtkSmartPointer int1Array = vtkSmartPointer::New(); int1Array->InsertNextValue(1); int1Array->InsertNextValue(0); int1Array->InsertNextValue(0); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(int1Array); vtkSmartPointer double1Array = vtkSmartPointer::New(); double1Array->InsertNextValue(center_1[0]); double1Array->InsertNextValue(center_1[1]); double1Array->InsertNextValue(center_1[2]); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(double1Array); double center_2[3] = {4.0f, 4.0f, 4.0f}; double normal_2[3] = {1.0f, 0.0f, 0.0f}; vtkSmartPointer p_source_2 = vtkSmartPointer::New(); p_source_2->SetNumberOfSides(80); p_source_2->SetCenter(center_2); p_source_2->SetRadius(4); p_source_2->SetNormal(normal_2); p_source_2->Update(); vtkPolyData *poly_2 = p_source_2->GetOutput(); mitk::Surface::Pointer surf_2 = mitk::Surface::New(); surf_2->SetVtkPolyData(poly_2); vtkSmartPointer int2Array = vtkSmartPointer::New(); int2Array->InsertNextValue(1); int2Array->InsertNextValue(0); int2Array->InsertNextValue(0); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(int2Array); vtkSmartPointer double2Array = vtkSmartPointer::New(); double2Array->InsertNextValue(center_2[0]); double2Array->InsertNextValue(center_2[1]); double2Array->InsertNextValue(center_2[2]); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(double2Array); std::vector surfaces; surfaces.push_back(surf_1); surfaces.push_back(surf_2); const mitk::PlaneGeometry * planeGeometry1 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry2 = GetPlaneGeometry(); std::vector planeGeometries; planeGeometries.push_back(planeGeometry1); planeGeometries.push_back(planeGeometry2); m_Controller->AddNewContours(surfaces, planeGeometries, true); // // Add contours CPPUNIT_ASSERT_MESSAGE("Wrong number of contours!", m_Controller->GetNumberOfContours() == 2); mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo3; contourInfo3.Contour = surf_1->Clone(); contourInfo3.ContourNormal = normal_1; contourInfo3.ContourPoint = center_1; // Shift the new contour so that it is different - contourInfo3.ContourPoint += 0.5; + contourInfo3.ContourPoint += mitk::Vector(0.5); bool success = m_Controller->RemoveContour(contourInfo3); CPPUNIT_ASSERT_MESSAGE("Remove failed - contour was unintentionally removed!", (m_Controller->GetNumberOfContours() == 2) && !success); mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo2; contourInfo2.ContourNormal = normal_2; contourInfo2.ContourPoint = center_2; contourInfo2.Contour = surf_2; success = m_Controller->RemoveContour(contourInfo2); CPPUNIT_ASSERT_MESSAGE("Remove failed - contour was not removed!", (m_Controller->GetNumberOfContours() == 1) && success); // // Let's see if the other contour No. 1 is still there - contourInfo3.ContourPoint -= 0.5; + contourInfo3.ContourPoint -= mitk::Vector(0.5); const mitk::Surface *remainingContour = m_Controller->GetContour(contourInfo3); CPPUNIT_ASSERT_MESSAGE( "Remove failed - contour was accidentally removed!", (m_Controller->GetNumberOfContours() == 1) && mitk::Equal(*(surf_1->GetVtkPolyData()), *(remainingContour->GetVtkPolyData()), 0.000001, true) && success); } bool AssertImagesEqual4D(mitk::LabelSetImage *img1, mitk::LabelSetImage *img2) { auto selector1 = mitk::ImageTimeSelector::New(); selector1->SetInput(img1); selector1->SetChannelNr(0); auto selector2 = mitk::ImageTimeSelector::New(); selector2->SetInput(img2); selector2->SetChannelNr(0); int numTs1 = img1->GetTimeSteps(); int numTs2 = img2->GetTimeSteps(); if (numTs1 != numTs2) { return false; } /*mitk::ImagePixelWriteAccessor accessor( img1 ); itk::Index<4> ind; ind[0] = 5; ind[1] = 5; ind[2] = 5; ind[3] = 2; accessor.SetPixelByIndex( ind, 7 );*/ for (int ts = 0; ts < numTs1; ++ts) { selector1->SetTimeNr(ts); selector2->SetTimeNr(ts); selector1->Update(); selector2->Update(); mitk::Image::Pointer imgSel1 = selector1->GetOutput(); mitk::Image::Pointer imgSel2 = selector2->GetOutput(); MITK_ASSERT_EQUAL(imgSel1, imgSel2, "Segmentation images are not equal"); } return true; } void TestSetCurrentInterpolationSession4D() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10, 5}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage4D(dimensions1); // mitk::Image * segmentationImage_1 = dynamic_cast(segmentation_1.GetPointer()); unsigned int dimensions2[] = {20, 10, 30, 4}; mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage4D(dimensions2); // Test 1 m_Controller->SetCurrentInterpolationSession(segmentation_1); auto currentSegmentation = dynamic_cast(m_Controller->GetCurrentSegmentation().GetPointer()); AssertImagesEqual4D(currentSegmentation, segmentation_1->Clone()); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_1.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); // Test 2 m_Controller->SetCurrentInterpolationSession(segmentation_2); // MITK_ASSERT_EQUAL(m_Controller->GetCurrentSegmentation(), segmentation_2->Clone(), "Segmentation images are not // equal"); currentSegmentation = dynamic_cast(m_Controller->GetCurrentSegmentation().GetPointer()); // AssertImagesEqual4D(currentSegmentation, segmentation_2->Clone()); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", currentSegmentation == segmentation_2.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test 3 m_Controller->SetCurrentInterpolationSession(segmentation_1); // MITK_ASSERT_EQUAL(m_Controller->GetCurrentSegmentation(), segmentation_1->Clone(), "Segmentation images are not // equal"); currentSegmentation = dynamic_cast(m_Controller->GetCurrentSegmentation().GetPointer()); AssertImagesEqual4D(currentSegmentation, segmentation_1->Clone()); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_1.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test 4 m_Controller->SetCurrentInterpolationSession(segmentation_1); // MITK_ASSERT_EQUAL(m_Controller->GetCurrentSegmentation(), segmentation_1->Clone(), "Segmentation images are not // equal"); currentSegmentation = dynamic_cast(m_Controller->GetCurrentSegmentation().GetPointer()); AssertImagesEqual4D(currentSegmentation, segmentation_1->Clone()); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_1.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test 5 m_Controller->SetCurrentInterpolationSession(nullptr); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().IsNull()); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); } void TestReplaceInterpolationSession4D() { // Create segmentation image unsigned int dimensions1[] = {10, 10, 10, 5}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage4D(dimensions1); m_Controller->SetCurrentInterpolationSession(segmentation_1); m_Controller->SetCurrentTimePoint(0); // Create some contours double center_1[3] = {1.25f, 3.43f, 4.44f}; double normal_1[3] = {0.25f, 1.76f, 0.93f}; vtkSmartPointer p_source = vtkSmartPointer::New(); p_source->SetNumberOfSides(20); p_source->SetCenter(center_1); p_source->SetRadius(4); p_source->SetNormal(normal_1); p_source->Update(); vtkPolyData *poly_1 = p_source->GetOutput(); mitk::Surface::Pointer surf_1 = mitk::Surface::New(); surf_1->SetVtkPolyData(poly_1); vtkSmartPointer int1Array = vtkSmartPointer::New(); int1Array->InsertNextValue(1); int1Array->InsertNextValue(0); int1Array->InsertNextValue(0); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(int1Array); vtkSmartPointer double1Array = vtkSmartPointer::New(); double1Array->InsertNextValue(center_1[0]); double1Array->InsertNextValue(center_1[1]); double1Array->InsertNextValue(center_1[2]); surf_1->GetVtkPolyData()->GetFieldData()->AddArray(double1Array); double center_2[3] = {4.0f, 4.0f, 4.0f}; double normal_2[3] = {1.0f, 0.0f, 0.0f}; vtkSmartPointer p_source_2 = vtkSmartPointer::New(); p_source_2->SetNumberOfSides(80); p_source_2->SetCenter(center_2); p_source_2->SetRadius(4); p_source_2->SetNormal(normal_2); p_source_2->Update(); vtkPolyData *poly_2 = p_source_2->GetOutput(); mitk::Surface::Pointer surf_2 = mitk::Surface::New(); surf_2->SetVtkPolyData(poly_2); vtkSmartPointer int2Array = vtkSmartPointer::New(); int2Array->InsertNextValue(1); int2Array->InsertNextValue(0); int2Array->InsertNextValue(0); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(int2Array); vtkSmartPointer double2Array = vtkSmartPointer::New(); double2Array->InsertNextValue(center_2[0]); double2Array->InsertNextValue(center_2[1]); double2Array->InsertNextValue(center_2[2]); surf_2->GetVtkPolyData()->GetFieldData()->AddArray(double2Array); std::vector surfaces; surfaces.push_back(surf_1); surfaces.push_back(surf_2); const mitk::PlaneGeometry * planeGeometry1 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry2 = GetPlaneGeometry(); std::vector planeGeometries; planeGeometries.push_back(planeGeometry1); planeGeometries.push_back(planeGeometry2); // Add contours m_Controller->AddNewContours(surfaces, planeGeometries, true); // Add contours for another timestep m_Controller->SetCurrentTimePoint(2); double center_3[3] = {1.3f, 3.5f, 4.6f}; double normal_3[3] = {0.20f, 1.6f, 0.8f}; vtkSmartPointer p_source_3 = vtkSmartPointer::New(); p_source_3->SetNumberOfSides(20); p_source_3->SetCenter(center_3); p_source_3->SetRadius(4); p_source_3->SetNormal(normal_3); p_source_3->Update(); vtkPolyData *poly_3 = p_source_3->GetOutput(); mitk::Surface::Pointer surf_3 = mitk::Surface::New(); surf_3->SetVtkPolyData(poly_3); vtkSmartPointer int3Array = vtkSmartPointer::New(); int3Array->InsertNextValue(1); int3Array->InsertNextValue(0); int3Array->InsertNextValue(2); surf_3->GetVtkPolyData()->GetFieldData()->AddArray(int3Array); vtkSmartPointer double3Array = vtkSmartPointer::New(); double3Array->InsertNextValue(center_3[0]); double3Array->InsertNextValue(center_3[1]); double3Array->InsertNextValue(center_3[2]); surf_3->GetVtkPolyData()->GetFieldData()->AddArray(double3Array); double center_4[3] = {1.32f, 3.53f, 4.8f}; double normal_4[3] = {0.22f, 1.5f, 0.85f}; vtkSmartPointer p_source_4 = vtkSmartPointer::New(); p_source_4->SetNumberOfSides(20); p_source_4->SetCenter(center_4); p_source_4->SetRadius(4); p_source_4->SetNormal(normal_4); p_source_4->Update(); vtkPolyData *poly_4 = p_source_4->GetOutput(); mitk::Surface::Pointer surf_4 = mitk::Surface::New(); surf_4->SetVtkPolyData(poly_4); vtkSmartPointer int4Array = vtkSmartPointer::New(); int4Array->InsertNextValue(1); int4Array->InsertNextValue(0); int4Array->InsertNextValue(2); surf_4->GetVtkPolyData()->GetFieldData()->AddArray(int4Array); vtkSmartPointer double4Array = vtkSmartPointer::New(); double4Array->InsertNextValue(center_4[0]); double4Array->InsertNextValue(center_4[1]); double4Array->InsertNextValue(center_4[2]); surf_4->GetVtkPolyData()->GetFieldData()->AddArray(double4Array); std::vector surfaces2; surfaces2.push_back(surf_3); surfaces2.push_back(surf_4); const mitk::PlaneGeometry * planeGeometry3 = GetPlaneGeometry(); const mitk::PlaneGeometry * planeGeometry4 = GetPlaneGeometry(); std::vector planeGeometries2; planeGeometries2.push_back(planeGeometry3); planeGeometries2.push_back(planeGeometry4); // Add contours m_Controller->AddNewContours(surfaces2, planeGeometries2, true); m_Controller->SetCurrentTimePoint(0); // Check if all contours are there mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo1; contourInfo1.ContourNormal = normal_1; contourInfo1.ContourPoint = center_1; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo2; contourInfo2.ContourNormal = normal_2; contourInfo2.ContourPoint = center_2; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo3; mitk::ScalarType n[3]; vtkPolygon::ComputeNormal(surf_3->GetVtkPolyData()->GetPoints(), n); contourInfo3.ContourNormal = n; contourInfo3.ContourPoint = center_3; mitk::SurfaceInterpolationController::ContourPositionInformation contourInfo4; // mitk::ScalarType n[3]; vtkPolygon::ComputeNormal(surf_4->GetVtkPolyData()->GetPoints(), n); contourInfo4.ContourNormal = n; contourInfo4.ContourPoint = center_4; const mitk::Surface *contour_1 = m_Controller->GetContour(contourInfo1); const mitk::Surface *contour_2 = m_Controller->GetContour(contourInfo2); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_1->GetVtkPolyData()), *(contour_1->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_2->GetVtkPolyData()), *(contour_2->GetVtkPolyData()), 0.000001, true)); m_Controller->SetCurrentTimePoint(2); const mitk::Surface *contour_3 = m_Controller->GetContour(contourInfo3); const mitk::Surface *contour_4 = m_Controller->GetContour(contourInfo4); CPPUNIT_ASSERT_MESSAGE("Wrong number of contours!", m_Controller->GetNumberOfContours() == 2); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_3->GetVtkPolyData()), *(contour_3->GetVtkPolyData()), 0.000001, true)); CPPUNIT_ASSERT_MESSAGE("Contours not equal!", mitk::Equal(*(surf_4->GetVtkPolyData()), *(contour_4->GetVtkPolyData()), 0.000001, true)); mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage4D(dimensions1); bool success = m_Controller->ReplaceInterpolationSession(segmentation_1, segmentation_2); CPPUNIT_ASSERT_MESSAGE("Replace session failed!", success == true); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_2.GetPointer()); } void TestRemoveAllInterpolationSessions4D() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10, 4}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage4D(dimensions1); unsigned int dimensions2[] = {20, 10, 30, 5}; mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage4D(dimensions2); // Test 1 m_Controller->SetCurrentInterpolationSession(segmentation_1); m_Controller->SetCurrentInterpolationSession(segmentation_2); m_Controller->RemoveAllInterpolationSessions(); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 0", m_Controller->GetNumberOfInterpolationSessions() == 0); } void TestRemoveInterpolationSession4D() { // Create image for testing unsigned int dimensions1[] = {10, 10, 10, 3}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage4D(dimensions1); unsigned int dimensions2[] = {20, 10, 30, 6}; mitk::LabelSetImage::Pointer segmentation_2 = createLabelSetImage4D(dimensions2); // Test 1 m_Controller->SetCurrentInterpolationSession(segmentation_1); m_Controller->SetCurrentInterpolationSession(segmentation_2); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test current segmentation should not be null if another one was removed m_Controller->RemoveInterpolationSession(segmentation_1); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); CPPUNIT_ASSERT_MESSAGE("Segmentation images are not equal", m_Controller->GetCurrentSegmentation().GetPointer() == segmentation_2.GetPointer()); CPPUNIT_ASSERT_MESSAGE("Current segmentation is null after another one was removed", m_Controller->GetCurrentSegmentation().IsNotNull()); m_Controller->SetCurrentInterpolationSession(segmentation_1); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 2", m_Controller->GetNumberOfInterpolationSessions() == 2); // Test current segmentation should not be null if another one was removed m_Controller->RemoveInterpolationSession(segmentation_1); CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 1", m_Controller->GetNumberOfInterpolationSessions() == 1); CPPUNIT_ASSERT_MESSAGE("Current segmentation is not null after session was removed", m_Controller->GetCurrentSegmentation().IsNull()); } void TestOnSegmentationDeleted4D() { { // Create image for testing unsigned int dimensions1[] = {10, 10, 10, 7}; mitk::LabelSetImage::Pointer segmentation_1 = createLabelSetImage4D(dimensions1); m_Controller->SetCurrentInterpolationSession(segmentation_1); m_Controller->SetCurrentTimePoint(3); } CPPUNIT_ASSERT_MESSAGE("Number of interpolation session not 0", m_Controller->GetNumberOfInterpolationSessions() == 0); } }; MITK_TEST_SUITE_REGISTRATION(mitkSurfaceInterpolationController) diff --git a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp index 87de1560a1..16952fd4a7 100644 --- a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp @@ -1,617 +1,617 @@ /*============================================================================ 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 "mitkCreateDistanceImageFromSurfaceFilter.h" #include "mitkImageCast.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkDoubleArray.h" #include "vtkPolyData.h" #include "vtkSmartPointer.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkNeighborhoodIterator.h" #include void mitk::CreateDistanceImageFromSurfaceFilter::CreateEmptyDistanceImage() { // Determine the bounds of the input points in index- and world-coordinates DistanceImageType::PointType minPointInWorldCoordinates, maxPointInWorldCoordinates; DistanceImageType::IndexType minPointInIndexCoordinates, maxPointInIndexCoordinates; DetermineBounds( minPointInWorldCoordinates, maxPointInWorldCoordinates, minPointInIndexCoordinates, maxPointInIndexCoordinates); // Calculate the extent of the region that contains all given points in MM. // To do this, we take the difference between the maximal and minimal // index-coordinates (must not be less than 1) and multiply it with the // spacing of the reference-image. Vector3D extentMM; for (unsigned int dim = 0; dim < 3; ++dim) { extentMM[dim] = (std::abs(maxPointInIndexCoordinates[dim] - minPointInIndexCoordinates[dim])) * m_ReferenceImage->GetSpacing()[dim]; } /* * Now create an empty distance image. The created image will always have the same number of pixels, independent from * the original image (e.g. always consists of 500000 pixels) and will have an isotropic spacing. * The spacing is calculated like the following: * The image's volume = 500000 Pixels = extentX*spacing*extentY*spacing*extentZ*spacing * So the spacing is: spacing = ( extentX*extentY*extentZ / 500000 )^(1/3) */ double basis = (extentMM[0] * extentMM[1] * extentMM[2]) / m_DistanceImageVolume; double exponent = 1.0 / 3.0; m_DistanceImageSpacing = pow(basis, exponent); // calculate the number of pixels of the distance image for each direction unsigned int numberOfXPixel = extentMM[0] / m_DistanceImageSpacing; unsigned int numberOfYPixel = extentMM[1] / m_DistanceImageSpacing; unsigned int numberOfZPixel = extentMM[2] / m_DistanceImageSpacing; // We increase the sizeOfRegion by 4 as we decrease the origin by 2 later. // This expansion of the region is necessary to achieve a complete // interpolation. DistanceImageType::SizeType sizeOfRegion; sizeOfRegion[0] = numberOfXPixel + 8; sizeOfRegion[1] = numberOfYPixel + 8; sizeOfRegion[2] = numberOfZPixel + 8; // The region starts at index 0,0,0 DistanceImageType::IndexType initialOriginAsIndex; initialOriginAsIndex.Fill(0); DistanceImageType::PointType originAsWorld = minPointInWorldCoordinates; DistanceImageType::RegionType lpRegion; lpRegion.SetSize(sizeOfRegion); lpRegion.SetIndex(initialOriginAsIndex); // We initialize the itk::Image with // * origin and direction to have it correctly placed and rotated in the world // * the largest possible region to set the extent to be calculated // * the isotropic spacing that we have calculated above m_DistanceImageITK = DistanceImageType::New(); m_DistanceImageITK->SetOrigin(originAsWorld); m_DistanceImageITK->SetDirection(m_ReferenceImage->GetDirection()); m_DistanceImageITK->SetRegions(lpRegion); - m_DistanceImageITK->SetSpacing(m_DistanceImageSpacing); + m_DistanceImageITK->SetSpacing(itk::Vector(m_DistanceImageSpacing)); m_DistanceImageITK->Allocate(); // First of all the image is initialized with the value 10*m_DistanceImageSpacing for each pixel m_DistanceImageDefaultBufferValue = 10 * m_DistanceImageSpacing; m_DistanceImageITK->FillBuffer(m_DistanceImageDefaultBufferValue); // Now we move the origin of the distanceImage 2 index-Coordinates // in all directions DistanceImageType::IndexType originAsIndex; m_DistanceImageITK->TransformPhysicalPointToIndex(originAsWorld, originAsIndex); originAsIndex[0] -= 2; originAsIndex[1] -= 2; originAsIndex[2] -= 2; m_DistanceImageITK->TransformIndexToPhysicalPoint(originAsIndex, originAsWorld); m_DistanceImageITK->SetOrigin(originAsWorld); } mitk::CreateDistanceImageFromSurfaceFilter::CreateDistanceImageFromSurfaceFilter() : m_DistanceImageSpacing(0.0), m_DistanceImageDefaultBufferValue(0.0) { m_DistanceImageVolume = 50000; this->m_UseProgressBar = false; this->m_ProgressStepSize = 5; mitk::Image::Pointer output = mitk::Image::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::CreateDistanceImageFromSurfaceFilter::~CreateDistanceImageFromSurfaceFilter() { } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateData() { this->PreprocessContourPoints(); this->CreateEmptyDistanceImage(); // First of all we have to build the equation-system from the existing contour-edge-points this->CreateSolutionMatrixAndFunctionValues(); if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(1); m_Weights = m_SolutionMatrix.partialPivLu().solve(m_FunctionValues); if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(2); // The last step is to create the distance map with the interpolated distance function this->FillDistanceImage(); if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(2); m_Centers.clear(); m_Normals.clear(); } void mitk::CreateDistanceImageFromSurfaceFilter::PreprocessContourPoints() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); if (numberOfInputs == 0) { MITK_ERROR << "mitk::CreateDistanceImageFromSurfaceFilter: No input available. Please set an input!" << std::endl; itkExceptionMacro("mitk::CreateDistanceImageFromSurfaceFilter: No input available. Please set an input!"); return; } // First of all we have to extract the nomals and the surface points. // Duplicated points can be eliminated vtkSmartPointer polyData; vtkSmartPointer currentCellNormals; vtkSmartPointer existingPolys; vtkSmartPointer existingPoints; double p[3]; PointType currentPoint; PointType normal; for (unsigned int i = 0; i < numberOfInputs; i++) { auto currentSurface = this->GetInput(i); polyData = currentSurface->GetVtkPolyData(); if (polyData->GetNumberOfPolys() == 0) { MITK_INFO << "mitk::CreateDistanceImageFromSurfaceFilter: No input-polygons available. Please be sure the input " "surface consists of polygons!" << std::endl; } currentCellNormals = vtkDoubleArray::SafeDownCast(polyData->GetCellData()->GetNormals()); existingPolys = polyData->GetPolys(); existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); const vtkIdType *cell(nullptr); vtkIdType cellSize(0); for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { for (vtkIdType j = 0; j < cellSize; j++) { existingPoints->GetPoint(cell[j], p); currentPoint.copy_in(p); int count = std::count(m_Centers.begin(), m_Centers.end(), currentPoint); if (count == 0) { double currentNormal[3]; currentCellNormals->GetTuple(cell[j], currentNormal); normal.copy_in(currentNormal); m_Normals.push_back(normal); m_Centers.push_back(currentPoint); } } // end for all points } // end for all cells } // end for all outputs } void mitk::CreateDistanceImageFromSurfaceFilter::CreateSolutionMatrixAndFunctionValues() { // For we can now calculate the exact size of the centers we initialize the data structures unsigned int numberOfCenters = m_Centers.size(); m_Centers.reserve(numberOfCenters * 3); m_FunctionValues.resize(numberOfCenters * 3); m_FunctionValues.fill(0); PointType currentPoint; PointType normal; // Create inner points for (unsigned int i = 0; i < numberOfCenters; i++) { currentPoint = m_Centers.at(i); normal = m_Normals.at(i); currentPoint[0] = currentPoint[0] - normal[0] * m_DistanceImageSpacing; currentPoint[1] = currentPoint[1] - normal[1] * m_DistanceImageSpacing; currentPoint[2] = currentPoint[2] - normal[2] * m_DistanceImageSpacing; m_Centers.push_back(currentPoint); m_FunctionValues[numberOfCenters + i] = -m_DistanceImageSpacing; } // Create outer points for (unsigned int i = 0; i < numberOfCenters; i++) { currentPoint = m_Centers.at(i); normal = m_Normals.at(i); currentPoint[0] = currentPoint[0] + normal[0] * m_DistanceImageSpacing; currentPoint[1] = currentPoint[1] + normal[1] * m_DistanceImageSpacing; currentPoint[2] = currentPoint[2] + normal[2] * m_DistanceImageSpacing; m_Centers.push_back(currentPoint); m_FunctionValues[numberOfCenters * 2 + i] = m_DistanceImageSpacing; } // Now we have created all centers and all function values. Next step is to create the solution matrix numberOfCenters = m_Centers.size(); m_SolutionMatrix.resize(numberOfCenters, numberOfCenters); m_Weights.resize(numberOfCenters); PointType p1; PointType p2; double norm; for (unsigned int i = 0; i < numberOfCenters; i++) { for (unsigned int j = 0; j < numberOfCenters; j++) { // Calculate the RBF value. Currently using Phi(r) = r with r is the euclidian distance between two points p1 = m_Centers.at(i); p2 = m_Centers.at(j); p1 = p1 - p2; norm = p1.two_norm(); m_SolutionMatrix(i, j) = norm; } } } void mitk::CreateDistanceImageFromSurfaceFilter::FillDistanceImage() { /* * Now we must calculate the distance for each pixel. But instead of calculating the distance value * for all of the image's pixels we proceed similar to the region growing algorithm: * * 1. Take the first pixel from the narrowband_point_list and calculate the distance for each neighbor (6er) * 2. If the current index's distance value is below a certain threshold push it into the list * 3. Next iteration take the next index from the list and originAsIndex with 1. again * * This is done until the narrowband_point_list is empty. */ typedef itk::ImageRegionIteratorWithIndex ImageIterator; typedef itk::NeighborhoodIterator NeighborhoodImageIterator; std::queue narrowbandPoints; PointType currentPoint = m_Centers.at(0); double distance = this->CalculateDistanceValue(currentPoint); // create itk::Point from vnl_vector DistanceImageType::PointType currentPointAsPoint; currentPointAsPoint[0] = currentPoint[0]; currentPointAsPoint[1] = currentPoint[1]; currentPointAsPoint[2] = currentPoint[2]; // Transform the input point in world-coordinates to index-coordinates DistanceImageType::IndexType currentIndex; m_DistanceImageITK->TransformPhysicalPointToIndex(currentPointAsPoint, currentIndex); assert( m_DistanceImageITK->GetLargestPossibleRegion().IsInside(currentIndex)); // we are quite certain this should hold narrowbandPoints.push(currentIndex); m_DistanceImageITK->SetPixel(currentIndex, distance); NeighborhoodImageIterator::RadiusType radius; radius.Fill(1); NeighborhoodImageIterator nIt(radius, m_DistanceImageITK, m_DistanceImageITK->GetLargestPossibleRegion()); unsigned int relativeNbIdx[] = {4, 10, 12, 14, 16, 22}; bool isInBounds = false; while (!narrowbandPoints.empty()) { nIt.SetLocation(narrowbandPoints.front()); narrowbandPoints.pop(); unsigned int *relativeNb = &relativeNbIdx[0]; for (int i = 0; i < 6; i++) { nIt.GetPixel(*relativeNb, isInBounds); if (isInBounds && nIt.GetPixel(*relativeNb) == m_DistanceImageDefaultBufferValue) { currentIndex = nIt.GetIndex(*relativeNb); // Transform the currently checked point from index-coordinates to // world-coordinates m_DistanceImageITK->TransformIndexToPhysicalPoint(currentIndex, currentPointAsPoint); // create a vnl_vector currentPoint[0] = currentPointAsPoint[0]; currentPoint[1] = currentPointAsPoint[1]; currentPoint[2] = currentPointAsPoint[2]; // and check the distance distance = this->CalculateDistanceValue(currentPoint); if (std::fabs(distance) <= m_DistanceImageSpacing * 2) { nIt.SetPixel(*relativeNb, distance); narrowbandPoints.push(currentIndex); } } relativeNb++; } } ImageIterator imgRegionIterator(m_DistanceImageITK, m_DistanceImageITK->GetLargestPossibleRegion()); imgRegionIterator.GoToBegin(); double prevPixelVal = 1; DistanceImageType::IndexType _size; _size.Fill(-1); _size += m_DistanceImageITK->GetLargestPossibleRegion().GetSize(); // Set every pixel inside the surface to -m_DistanceImageDefaultBufferValue except the edge point (so that the // received surface is closed) while (!imgRegionIterator.IsAtEnd()) { if (imgRegionIterator.Get() == m_DistanceImageDefaultBufferValue && prevPixelVal < 0) { while (imgRegionIterator.Get() == m_DistanceImageDefaultBufferValue) { if (imgRegionIterator.GetIndex()[0] == _size[0] || imgRegionIterator.GetIndex()[1] == _size[1] || imgRegionIterator.GetIndex()[2] == _size[2] || imgRegionIterator.GetIndex()[0] == 0U || imgRegionIterator.GetIndex()[1] == 0U || imgRegionIterator.GetIndex()[2] == 0U) { imgRegionIterator.Set(m_DistanceImageDefaultBufferValue); prevPixelVal = m_DistanceImageDefaultBufferValue; ++imgRegionIterator; break; } else { imgRegionIterator.Set((-1) * m_DistanceImageDefaultBufferValue); ++imgRegionIterator; prevPixelVal = (-1) * m_DistanceImageDefaultBufferValue; } } } else if (imgRegionIterator.GetIndex()[0] == _size[0] || imgRegionIterator.GetIndex()[1] == _size[1] || imgRegionIterator.GetIndex()[2] == _size[2] || imgRegionIterator.GetIndex()[0] == 0U || imgRegionIterator.GetIndex()[1] == 0U || imgRegionIterator.GetIndex()[2] == 0U) { imgRegionIterator.Set(m_DistanceImageDefaultBufferValue); prevPixelVal = m_DistanceImageDefaultBufferValue; ++imgRegionIterator; } else { prevPixelVal = imgRegionIterator.Get(); ++imgRegionIterator; } } Image::Pointer resultImage = this->GetOutput(); // Cast the created distance-Image from itk::Image to the mitk::Image // that is our output. CastToMitkImage(m_DistanceImageITK, resultImage); } double mitk::CreateDistanceImageFromSurfaceFilter::CalculateDistanceValue(PointType p) { double distanceValue(0); PointType p1; PointType p2; double norm; CenterList::iterator centerIter; unsigned int count(0); for (centerIter = m_Centers.begin(); centerIter != m_Centers.end(); centerIter++) { p1 = *centerIter; p2 = p - p1; norm = p2.two_norm(); distanceValue = distanceValue + (norm * m_Weights[count]); ++count; } return distanceValue; } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateOutputInformation() { } void mitk::CreateDistanceImageFromSurfaceFilter::PrintEquationSystem() { std::stringstream out; out << "Nummber of rows: " << m_SolutionMatrix.rows() << " ****** Number of columns: " << m_SolutionMatrix.cols() << endl; out << "[ "; for (int i = 0; i < m_SolutionMatrix.rows(); i++) { for (int j = 0; j < m_SolutionMatrix.cols(); j++) { out << m_SolutionMatrix(i, j) << " "; } out << ";" << endl; } out << " ]\n\n\n"; for (unsigned int i = 0; i < m_Centers.size(); i++) { out << m_Centers.at(i) << ";" << endl; } std::cout << "Equation system: \n\n\n" << out.str(); } void mitk::CreateDistanceImageFromSurfaceFilter::SetInput(const mitk::Surface *surface) { this->SetInput(0, surface); } void mitk::CreateDistanceImageFromSurfaceFilter::SetInput(unsigned int idx, const mitk::Surface *surface) { if (this->GetInput(idx) != surface) { this->SetNthInput(idx, const_cast(surface)); this->Modified(); } } const mitk::Surface *mitk::CreateDistanceImageFromSurfaceFilter::GetInput() { if (this->GetNumberOfIndexedInputs() < 1) return nullptr; return static_cast(this->ProcessObject::GetInput(0)); } const mitk::Surface *mitk::CreateDistanceImageFromSurfaceFilter::GetInput(unsigned int idx) { if (this->GetNumberOfIndexedInputs() < 1) return nullptr; return static_cast(this->ProcessObject::GetInput(idx)); } void mitk::CreateDistanceImageFromSurfaceFilter::RemoveInputs(mitk::Surface *input) { DataObjectPointerArraySizeType nb = this->GetNumberOfIndexedInputs(); for (DataObjectPointerArraySizeType i = 0; i < nb; i++) { if (this->GetInput(i) == input) { this->RemoveInput(i); return; } } } void mitk::CreateDistanceImageFromSurfaceFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(1); mitk::Image::Pointer output = mitk::Image::New(); this->SetNthOutput(0, output.GetPointer()); } void mitk::CreateDistanceImageFromSurfaceFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::CreateDistanceImageFromSurfaceFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } void mitk::CreateDistanceImageFromSurfaceFilter::SetReferenceImage(itk::ImageBase<3>::Pointer referenceImage) { m_ReferenceImage = referenceImage; } void mitk::CreateDistanceImageFromSurfaceFilter::DetermineBounds( DistanceImageType::PointType &minPointInWorldCoordinates, DistanceImageType::PointType &maxPointInWorldCoordinates, DistanceImageType::IndexType &minPointInIndexCoordinates, DistanceImageType::IndexType &maxPointInIndexCoordinates) { PointType firstCenter = m_Centers.at(0); DistanceImageType::PointType tmpPoint; tmpPoint[0] = firstCenter[0]; tmpPoint[1] = firstCenter[1]; tmpPoint[2] = firstCenter[2]; // transform the first point from world-coordinates to index-coordinates itk::ContinuousIndex tmpIndex; m_ReferenceImage->TransformPhysicalPointToContinuousIndex(tmpPoint, tmpIndex); // initialize the variables with this first point DistanceImageType::IndexValueType xmin = tmpIndex[0]; DistanceImageType::IndexValueType ymin = tmpIndex[1]; DistanceImageType::IndexValueType zmin = tmpIndex[2]; DistanceImageType::IndexValueType xmax = tmpIndex[0]; DistanceImageType::IndexValueType ymax = tmpIndex[1]; DistanceImageType::IndexValueType zmax = tmpIndex[2]; // iterate over the rest of the points auto centerIter = m_Centers.begin(); for (++centerIter; centerIter != m_Centers.end(); centerIter++) { tmpPoint[0] = (*centerIter)[0]; tmpPoint[1] = (*centerIter)[1]; tmpPoint[2] = (*centerIter)[2]; // transform each point from world-coordinates to index-coordinates m_ReferenceImage->TransformPhysicalPointToContinuousIndex(tmpPoint, tmpIndex); // and set the variables accordingly to find the minimum // and maximum in all directions in index-coordinates if (xmin > tmpIndex[0]) { xmin = tmpIndex[0]; } if (ymin > tmpIndex[1]) { ymin = tmpIndex[1]; } if (zmin > tmpIndex[2]) { zmin = tmpIndex[2]; } if (xmax < tmpIndex[0]) { xmax = tmpIndex[0]; } if (ymax < tmpIndex[1]) { ymax = tmpIndex[1]; } if (zmax < tmpIndex[2]) { zmax = tmpIndex[2]; } } // put the found coordinates into Index-Points minPointInIndexCoordinates[0] = xmin; minPointInIndexCoordinates[1] = ymin; minPointInIndexCoordinates[2] = zmin; maxPointInIndexCoordinates[0] = xmax; maxPointInIndexCoordinates[1] = ymax; maxPointInIndexCoordinates[2] = zmax; // and transform them into world-coordinates m_ReferenceImage->TransformIndexToPhysicalPoint(minPointInIndexCoordinates, minPointInWorldCoordinates); m_ReferenceImage->TransformIndexToPhysicalPoint(maxPointInIndexCoordinates, maxPointInWorldCoordinates); } diff --git a/Plugins/org.mitk.gui.qt.tubegraph/src/internal/QmitkTubeGraphLabelWidget.cpp b/Plugins/org.mitk.gui.qt.tubegraph/src/internal/QmitkTubeGraphLabelWidget.cpp index 05776510ec..ca5289bc97 100644 --- a/Plugins/org.mitk.gui.qt.tubegraph/src/internal/QmitkTubeGraphLabelWidget.cpp +++ b/Plugins/org.mitk.gui.qt.tubegraph/src/internal/QmitkTubeGraphLabelWidget.cpp @@ -1,109 +1,109 @@ /*============================================================================ 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 "QmitkTubeGraphLabelWidget.h" #include #include QmitkTubeGraphLabelWidget::QmitkTubeGraphLabelWidget(QWidget* parent) :QWidget(parent) { this->InitWidget(); } QmitkTubeGraphLabelWidget::~QmitkTubeGraphLabelWidget() { delete m_VisibilityCheckBox; delete m_LabelPushButton; delete m_ColoringPushButton; } void QmitkTubeGraphLabelWidget::InitWidget() { m_VisibilityCheckBox = new QCheckBox("", this); m_VisibilityCheckBox->setChecked(true); m_VisibilityCheckBox->setSizePolicy(QSizePolicy::Fixed, QSizePolicy::Fixed); m_LabelPushButton = new QPushButton("", this); //m_LabelPushButton->setCheckable(true); //m_LabelPushButton->setAutoExclusive(true); m_ColoringPushButton = new QPushButton(this); m_ColoringPushButton->setStyleSheet( QString("* { background-color: rgb(0,0,0)}")); m_ColoringPushButton->setSizePolicy(QSizePolicy::Fixed, QSizePolicy::Fixed); auto layout = new QHBoxLayout; layout->addWidget(m_VisibilityCheckBox); layout->addWidget(m_LabelPushButton); layout->addWidget(m_ColoringPushButton); this->setLayout(layout); connect(m_VisibilityCheckBox, SIGNAL(toggled(bool)), this, SLOT(OnVisibilityToggled(bool))); connect(m_LabelPushButton, SIGNAL(clicked()), this, SLOT(OnLabelButtonClicked())); connect(m_ColoringPushButton, SIGNAL(clicked()), this, SLOT(OnColoringButtonClicked())); //connect(m_ColoringPushButton, SIGNAL(clicked()), this, SIGNAL(SignalLabelColorChanged(this->GetLabelName()))); } void QmitkTubeGraphLabelWidget::SetLabelName(QString name) { m_LabelPushButton->setText(name); } QString QmitkTubeGraphLabelWidget::GetLabelName() { return m_LabelPushButton->text(); } void QmitkTubeGraphLabelWidget::SetLabelColor(mitk::Color* color) { m_ColoringPushButton->setStyleSheet( QString("* { background-color: rgb(%1,%2,%3)}").arg(color->GetRed()).arg(color->GetGreen()).arg(color->GetBlue()) ); } mitk::Color* QmitkTubeGraphLabelWidget::GetLabelColor() { //return the background color of the button const QColor color = m_ColoringPushButton->palette().color(QPalette::Button); - auto rgb = new mitk::Color(); - rgb[0] = color.red(); rgb[1] = color.green(); rgb[2] = color.blue(); + auto rgb = new mitk::Color(); + rgb->Set(color.red(), color.green(), color.blue()); return rgb; } void QmitkTubeGraphLabelWidget::OnVisibilityToggled(bool isVisible) { emit SignalLabelVisibilityToggled(isVisible, this->GetLabelName()); } void QmitkTubeGraphLabelWidget::OnLabelButtonClicked() { emit SignalLabelButtonClicked(this->GetLabelName()); } void QmitkTubeGraphLabelWidget::OnColoringButtonClicked() { QColor usersChoice = QColorDialog::getColor(); mitk::Color color; if (usersChoice.spec() != 0) { color[0] = usersChoice.red(); color[1] = usersChoice.green(); color[2] = usersChoice.blue(); this->SetLabelColor(&color); emit SignalLabelColorChanged(color, this->GetLabelName()); } }