diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp index 8b32f1f143..36c0894ac8 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp @@ -1,375 +1,361 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTrackingHandlerTensor.h" namespace mitk { TrackingHandlerTensor::TrackingHandlerTensor() : m_FaThreshold(0.1) , m_F(1.0) , m_G(0.0) , m_InterpolateTensors(true) , m_NumberOfInputs(0) { } TrackingHandlerTensor::~TrackingHandlerTensor() { } void TrackingHandlerTensor::InitForTracking() { MITK_INFO << "Initializing tensor tracker."; if (m_NeedsDataInit) { m_NumberOfInputs = m_TensorImages.size(); for (int i=0; iSetSpacing( m_TensorImages.at(0)->GetSpacing() ); pdImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); pdImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); pdImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); pdImage->Allocate(); m_PdImage.push_back(pdImage); ItkDoubleImgType::Pointer emaxImage = ItkDoubleImgType::New(); emaxImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); emaxImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); emaxImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); emaxImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); emaxImage->Allocate(); emaxImage->FillBuffer(1.0); m_EmaxImage.push_back(emaxImage); } bool useUserFaImage = true; if (m_FaImage.IsNull()) { m_FaImage = ItkFloatImgType::New(); m_FaImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); m_FaImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); m_FaImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); m_FaImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); m_FaImage->Allocate(); m_FaImage->FillBuffer(0.0); useUserFaImage = false; } typedef itk::DiffusionTensor3D TensorType; - TensorType::EigenValuesArrayType eigenvalues; - TensorType::EigenVectorsMatrixType eigenvectors; - -#ifdef WIN32 -#pragma omp parallel for -#else -#pragma omp parallel for collapse(3) -#endif + for (int x=0; x<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[0]; x++) for (int y=0; y<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[1]; y++) for (int z=0; z<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[2]; z++) { ItkTensorImageType::IndexType index; index[0] = x; index[1] = y; index[2] = z; for (int i=0; iGetPixel(index); - tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); - } - vnl_vector_fixed dir; + tensor = m_TensorImages.at(i)->GetPixel(index); + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude()>mitk::eps) dir.normalize(); else dir.fill(0.0); -#pragma omp critical - { - m_PdImage.at(i)->SetPixel(index, dir); - if (!useUserFaImage) - m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy()); - m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]); - } + m_PdImage.at(i)->SetPixel(index, dir); + if (!useUserFaImage) + m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy()); + m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]); } if (!useUserFaImage) - { -#pragma omp critical m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)/m_NumberOfInputs); - } } m_NeedsDataInit = false; } if (m_F+m_G>1.0) { float temp = m_F+m_G; m_F /= temp; m_G /= temp; } std::cout << "TrackingHandlerTensor - FA threshold: " << m_FaThreshold << std::endl; std::cout << "TrackingHandlerTensor - f: " << m_F << std::endl; std::cout << "TrackingHandlerTensor - g: " << m_G << std::endl; } vnl_vector_fixed TrackingHandlerTensor::GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magGetPixel(idx); if (out_dir.magnitude()>0.5) { image_num = i; oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } else { for (unsigned int i=0; i dir = m_PdImage.at(i)->GetPixel(idx); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { image_num = i; angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } vnl_vector_fixed TrackingHandlerTensor::GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor) { // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; m_FaImage->TransformPhysicalPointToIndex(itkP, idx); m_FaImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_FaImage->GetLargestPossibleRegion().IsInside(idx) ) return dir; int image_num = -1; if (!m_Interpolate) { dir = GetMatchingDirection(idx, oldDir, image_num); if (image_num>=0) tensor = m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx); } else { float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); dir = GetMatchingDirection(idx, oldDir, image_num) * interpWeights[0]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx) * interpWeights[0]; itk::Index<3> tmpIdx = idx; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[1]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[2]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[3]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[4]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[5]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[6]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[7]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[7]; } } return dir; } vnl_vector_fixed TrackingHandlerTensor::GetLargestEigenvector(TensorType& tensor) { vnl_vector_fixed dir; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude() TrackingHandlerTensor::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); TensorType tensor; tensor.Fill(0); try { itk::Index<3> index; m_TensorImages.at(0)->TransformPhysicalPointToIndex(pos, index); float fa = GetImageValue(pos, m_FaImage, m_Interpolate); if (fa oldDir = olddirs.back(); if (m_FlipX) oldDir[0] *= -1; if (m_FlipY) oldDir[1] *= -1; if (m_FlipZ) oldDir[2] *= -1; float old_mag = oldDir.magnitude(); if (!m_Interpolate && oldIndex==index) return oldDir; output_direction = GetDirection(pos, oldDir, tensor); float mag = output_direction.magnitude(); if (mag>=mitk::eps) { output_direction.normalize(); if (old_mag>0.5 && m_G>mitk::eps) // TEND tracking { output_direction[0] = m_F*output_direction[0] + (1-m_F)*( (1-m_G)*oldDir[0] + m_G*(tensor[0]*oldDir[0] + tensor[1]*oldDir[1] + tensor[2]*oldDir[2])); output_direction[1] = m_F*output_direction[1] + (1-m_F)*( (1-m_G)*oldDir[1] + m_G*(tensor[1]*oldDir[0] + tensor[3]*oldDir[1] + tensor[4]*oldDir[2])); output_direction[2] = m_F*output_direction[2] + (1-m_F)*( (1-m_G)*oldDir[2] + m_G*(tensor[2]*oldDir[0] + tensor[4]*oldDir[1] + tensor[5]*oldDir[2])); output_direction.normalize(); } float a = 1; if (old_mag>0.5) a = dot_product(output_direction, oldDir); if (a>=m_AngularThreshold) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); } catch(...) { } if (m_FlipX) output_direction[0] *= -1; if (m_FlipY) output_direction[1] *= -1; if (m_FlipZ) output_direction[2] *= -1; return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h index b0957d7c6c..cf5b2455df 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h @@ -1,88 +1,88 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingHandlerTensor #define _TrackingHandlerTensor #include "mitkTrackingDataHandler.h" #include #include namespace mitk { /** * \brief Enables streamline tracking on tensor images. Supports multi tensor tracking by adding multiple tensor images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerTensor : public TrackingDataHandler { public: TrackingHandlerTensor(); ~TrackingHandlerTensor(); typedef itk::DiffusionTensor3D TensorType; typedef itk::Image< TensorType, 3 > ItkTensorImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position void SetF(float f){ m_F = f; } void SetG(float g){ m_G = g; } void SetFaThreshold(float FaThreshold){ m_FaThreshold = FaThreshold; } - void AddTensorImage( ItkTensorImageType::Pointer img ){ m_TensorImages.push_back(img); DataModified(); } - void SetTensorImage( ItkTensorImageType::Pointer img ){ m_TensorImages.clear(); m_TensorImages.push_back(img); DataModified(); } + void AddTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.push_back(img); DataModified(); } + void SetTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.clear(); m_TensorImages.push_back(img); DataModified(); } void ClearTensorImages(){ m_TensorImages.clear(); DataModified(); } void SetFaImage( ItkFloatImgType::Pointer img ){ m_FaImage = img; DataModified(); } void SetInterpolateTensors( bool interpolateTensors ){ m_InterpolateTensors = interpolateTensors; } void SetMode( MODE m ) { if (m==MODE::DETERMINISTIC) m_Mode = m; else mitkThrow() << "Tensor tracker is only implemented for deterministic mode."; } ItkUcharImgType::SpacingType GetSpacing(){ return m_FaImage->GetSpacing(); } itk::Point GetOrigin(){ return m_FaImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection(){ return m_FaImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion(){ return m_FaImage->GetLargestPossibleRegion(); } protected: vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num); vnl_vector_fixed GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor); vnl_vector_fixed GetLargestEigenvector(TensorType& tensor); float m_FaThreshold; float m_F; float m_G; std::vector< ItkDoubleImgType::Pointer > m_EmaxImage; ///< Stores largest eigenvalues per voxel (one for each tensor) ItkFloatImgType::Pointer m_FaImage; ///< FA image used to determine streamline termination. std::vector< ItkPDImgType::Pointer > m_PdImage; ///< Stores principal direction of each tensor in each voxel. - std::vector< ItkTensorImageType::Pointer > m_TensorImages; ///< Input tensor images. For multi tensor tracking provide multiple tensor images. + std::vector< ItkTensorImageType::ConstPointer > m_TensorImages; ///< Input tensor images. For multi tensor tracking provide multiple tensor images. bool m_InterpolateTensors; ///< If false, then the peaks are interpolated. Otherwiese, The tensors are interpolated. int m_NumberOfInputs; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp index 344f3ee3a6..5495d086b1 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp @@ -1,429 +1,429 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class mitkStreamlineTractographyTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkStreamlineTractographyTestSuite); MITK_TEST(Test_Peak1); MITK_TEST(Test_Peak2); MITK_TEST(Test_Tensor1); MITK_TEST(Test_Tensor2); MITK_TEST(Test_Tensor3); MITK_TEST(Test_Odf1); MITK_TEST(Test_Odf2); MITK_TEST(Test_Odf3); MITK_TEST(Test_Odf4); MITK_TEST(Test_Odf5); MITK_TEST(Test_Odf6); CPPUNIT_TEST_SUITE_END(); typedef itk::VectorImage< short, 3> ItkDwiType; private: public: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ typedef itk::Image ItkUcharImgType; typedef itk::Image ItkFloatImgType; mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itk_qball_image; - mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itk_tensor_image; + mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itk_tensor_image; mitk::TrackingHandlerPeaks::PeakImgType::Pointer itk_peak_image; ItkUcharImgType::Pointer itk_seed_image; ItkUcharImgType::Pointer itk_mask_image; ItkFloatImgType::Pointer itk_gfa_image; float gfa_threshold; float odf_threshold; float peak_threshold; itk::StreamlineTrackingFilter::Pointer tracker; void setUp() override { omp_set_num_threads(1); gfa_threshold = 0.2; odf_threshold = 0.1; peak_threshold = 0.1; mitk::Image::Pointer qball_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/qball_image.qbi")); mitk::Image::Pointer tensor_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/tensor_image.dti")); mitk::Image::Pointer peak_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/qball_peak_image.nii.gz")); mitk::Image::Pointer seed_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/seed_image.nii.gz")); mitk::Image::Pointer mask_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/mask_image.nii.gz")); mitk::Image::Pointer gfa_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/gfa_image.nii.gz")); { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(peak_image); caster->Update(); itk_peak_image = caster->GetOutput(); } { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(tensor_image); caster->Update(); itk_tensor_image = caster->GetOutput(); } { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(qball_image); caster->Update(); itk_qball_image = caster->GetOutput(); } itk_gfa_image = ItkFloatImgType::New(); mitk::CastToItkImage(gfa_image, itk_gfa_image); itk_seed_image = ItkUcharImgType::New(); mitk::CastToItkImage(seed_image, itk_seed_image); itk_mask_image = ItkUcharImgType::New(); mitk::CastToItkImage(mask_image, itk_mask_image); } mitk::FiberBundle::Pointer LoadReferenceFib(std::string filename) { mitk::FiberBundle::Pointer fib = nullptr; if (itksys::SystemTools::FileExists(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename))) { mitk::BaseData::Pointer baseData = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename)).at(0); fib = dynamic_cast(baseData.GetPointer()); } return fib; } mitk::Image::Pointer LoadReferenceImage(std::string filename) { mitk::Image::Pointer img = nullptr; if (itksys::SystemTools::FileExists(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename))) { img = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename)); } return img; } void SetupTracker(mitk::TrackingDataHandler* handler) { tracker = itk::StreamlineTrackingFilter::New(); tracker->SetRandom(false); tracker->SetNumberOfSamples(0); tracker->SetAngularThreshold(-1); tracker->SetMaskImage(itk_mask_image); tracker->SetSeedImage(itk_seed_image); tracker->SetStoppingRegions(nullptr); tracker->SetSeedsPerVoxel(1); tracker->SetStepSize(0.5); tracker->SetSamplingDistance(0.25); tracker->SetUseStopVotes(true); tracker->SetOnlyForwardSamples(true); tracker->SetAposterioriCurvCheck(false); tracker->SetTissueImage(nullptr); tracker->SetSeedOnlyGm(false); tracker->SetControlGmEndings(false); tracker->SetMinTractLength(20); tracker->SetMaxNumTracts(-1); tracker->SetTrackingHandler(handler); tracker->SetUseOutputProbabilityMap(false); } void tearDown() override { } void CheckFibResult(std::string ref_file, mitk::FiberBundle::Pointer test_fib) { mitk::FiberBundle::Pointer ref = LoadReferenceFib(ref_file); if (ref.IsNull()) { mitk::IOUtil::Save(test_fib, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Reference file not found. Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } else { bool is_equal = ref->Equals(test_fib); if (!is_equal) { mitk::IOUtil::Save(test_fib, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Tractograms are not equal! Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } } } void CheckImageResult(std::string ref_file, mitk::Image::Pointer test_img) { mitk::Image::Pointer ref = LoadReferenceImage(ref_file); if (ref.IsNull()) { mitk::IOUtil::Save(test_img, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Reference file not found. Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } else { MITK_ASSERT_EQUAL(test_img, ref, "Images should be equal"); } } void Test_Peak1() { mitk::TrackingHandlerPeaks* handler = new mitk::TrackingHandlerPeaks(); handler->SetPeakImage(itk_peak_image); handler->SetPeakThreshold(peak_threshold); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Peak1.fib", outFib); delete handler; } void Test_Peak2() { mitk::TrackingHandlerPeaks* handler = new mitk::TrackingHandlerPeaks(); handler->SetPeakImage(itk_peak_image); handler->SetPeakThreshold(peak_threshold); handler->SetInterpolate(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Peak2.fib", outFib); delete handler; } void Test_Tensor1() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); handler->SetFaThreshold(gfa_threshold); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor1.fib", outFib); delete handler; } void Test_Tensor2() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); handler->SetFaThreshold(gfa_threshold); handler->SetInterpolate(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor2.fib", outFib); delete handler; } void Test_Tensor3() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); handler->SetFaThreshold(gfa_threshold); handler->SetInterpolate(false); handler->SetF(0); handler->SetG(1); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor3.fib", outFib); delete handler; } void Test_Odf1() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_qball_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf1.fib", outFib); delete handler; } void Test_Odf2() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_qball_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf2.fib", outFib); delete handler; } void Test_Odf3() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_qball_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); handler->SetInterpolate(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf3.fib", outFib); delete handler; } void Test_Odf4() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_qball_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); SetupTracker(handler); tracker->SetSeedsPerVoxel(3); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf4.fib", outFib); delete handler; } void Test_Odf5() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_qball_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); SetupTracker(handler); tracker->SetSeedsPerVoxel(3); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf5.fib", outFib); delete handler; } void Test_Odf6() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_qball_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); SetupTracker(handler); tracker->SetSeedsPerVoxel(10); tracker->SetUseOutputProbabilityMap(true); tracker->Update(); itk::StreamlineTrackingFilter::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); CheckImageResult("Test_Odf6.nrrd", img); delete handler; } }; MITK_TEST_SUITE_REGISTRATION(mitkStreamlineTractography) diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp index 99966c5927..0c13c89850 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp @@ -1,471 +1,471 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include using namespace std; const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false); parser.addArgument("algorithm", "a", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle/probability map", us::Any(), false); parser.addArgument("stop_mask", "", mitkCommandLineParser::String, "Stop image:", "streamlines entering the binary mask will stop immediately", us::Any()); parser.addArgument("tracking_mask", "", mitkCommandLineParser::String, "Mask image:", "restrict tractography with a binary mask image", us::Any()); parser.addArgument("seed_mask", "", mitkCommandLineParser::String, "Seed image:", "binary mask image defining seed voxels", us::Any()); parser.addArgument("tissue_image", "", mitkCommandLineParser::String, "Tissue type image:", "image with tissue type labels (WM=3, GM=1)", us::Any()); parser.addArgument("sharpen_odfs", "", mitkCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkCommandLineParser::Float, "ODF Cutoff:", "additional threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.1); parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size)"); parser.addArgument("min_tract_length", "", mitkCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_gm", "", mitkCommandLineParser::Bool, "Seed only GM:", "Seed only in gray matter (requires tissue type image --tissue_image)"); parser.addArgument("control_gm_endings", "", mitkCommandLineParser::Bool, "Control GM endings:", "Seed perpendicular to gray matter and enforce endings inside of the gray matter (requires tissue type image --tissue_image)"); parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded.", -1); parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.addArgument("output_prob_map", "", mitkCommandLineParser::Bool, "Output probability map:", "output probability map instead of tractogram"); parser.addArgument("no_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate:", "don't interpolate image values"); parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); //parser.addArgument("apply_image_rotation", "", mitkCommandLineParser::Bool, "Apply image rotation:", "applies image rotation to image peaks (only for 'Peaks' algorithm). This is necessary in some cases, e.g. if the peaks were obtained with MRtrix."); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "Compress output fibers using the given error threshold (in mm)"); parser.addArgument("additional_images", "", mitkCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional Features)", us::Any()); // parameters for random forest based tractography parser.addArgument("forest", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any()); parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); // parameters for tensor tractography parser.addArgument("tend_f", "", mitkCommandLineParser::Float, "Weight f", "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkCommandLineParser::Float, "Weight g", "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); string outFile = us::any_cast(parsedArgs["out"]); string algorithm = us::any_cast(parsedArgs["algorithm"]); bool sharpen_odfs = false; if (parsedArgs.count("sharpen_odfs")) sharpen_odfs = us::any_cast(parsedArgs["sharpen_odfs"]); bool interpolate = true; if (parsedArgs.count("no_interpolation")) interpolate = !us::any_cast(parsedArgs["no_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); bool seed_gm = false; if (parsedArgs.count("seed_gm")) seed_gm = us::any_cast(parsedArgs["seed_gm"]); bool control_gm_endings = false; if (parsedArgs.count("control_gm_endings")) control_gm_endings = us::any_cast(parsedArgs["control_gm_endings"]); bool use_stop_votes = false; if (parsedArgs.count("use_stop_votes")) use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); bool use_only_forward_samples = false; if (parsedArgs.count("use_only_forward_samples")) use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); bool output_prob_map = false; if (parsedArgs.count("output_prob_map")) output_prob_map = us::any_cast(parsedArgs["output_prob_map"]); bool flip_x = false; if (parsedArgs.count("flip_x")) flip_x = us::any_cast(parsedArgs["flip_x"]); bool flip_y = false; if (parsedArgs.count("flip_y")) flip_y = us::any_cast(parsedArgs["flip_y"]); bool flip_z = false; if (parsedArgs.count("flip_z")) flip_z = us::any_cast(parsedArgs["flip_z"]); bool apply_image_rotation = false; if (parsedArgs.count("apply_image_rotation")) apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); float min_tract_length = 20; if (parsedArgs.count("min_tract_length")) min_tract_length = us::any_cast(parsedArgs["min_tract_length"]); string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); string seedFile = ""; if (parsedArgs.count("seed_mask")) seedFile = us::any_cast(parsedArgs["seed_mask"]); string stopFile = ""; if (parsedArgs.count("stop_mask")) stopFile = us::any_cast(parsedArgs["stop_mask"]); string tissueFile = ""; if (parsedArgs.count("tissue_image")) tissueFile = us::any_cast(parsedArgs["tissue_image"]); float cutoff = 0.1; if (parsedArgs.count("cutoff")) cutoff = us::any_cast(parsedArgs["cutoff"]); float odf_cutoff = 0.1; if (parsedArgs.count("odf_cutoff")) odf_cutoff = us::any_cast(parsedArgs["odf_cutoff"]); float stepsize = -1; if (parsedArgs.count("step_size")) stepsize = us::any_cast(parsedArgs["step_size"]); float sampling_distance = -1; if (parsedArgs.count("sampling_distance")) sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); int num_samples = 0; if (parsedArgs.count("num_samples")) num_samples = us::any_cast(parsedArgs["num_samples"]); int seeds = 1; if (parsedArgs.count("seeds")) seeds = us::any_cast(parsedArgs["seeds"]); float tend_f = 1; if (parsedArgs.count("tend_f")) tend_f = us::any_cast(parsedArgs["tend_f"]); float tend_g = 0; if (parsedArgs.count("tend_g")) tend_g = us::any_cast(parsedArgs["tend_g"]); float angular_threshold = -1; if (parsedArgs.count("angular_threshold")) angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); unsigned int max_tracts = -1; if (parsedArgs.count("max_tracts")) max_tracts = us::any_cast(parsedArgs["max_tracts"]); // LOAD DATASETS mitkCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkUcharImgType; MITK_INFO << "loading input"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(mitk::IOUtil::Load(input_files.at(i))[0].GetPointer()); input_images.push_back(mitkImage); } ItkUcharImgType::Pointer mask; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(maskFile)[0].GetPointer()); mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); } ItkUcharImgType::Pointer seed; if (!seedFile.empty()) { MITK_INFO << "loading seed image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(seedFile)[0].GetPointer()); seed = ItkUcharImgType::New(); mitk::CastToItkImage(img, seed); } ItkUcharImgType::Pointer stop; if (!stopFile.empty()) { MITK_INFO << "loading stop image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(stopFile)[0].GetPointer()); stop = ItkUcharImgType::New(); mitk::CastToItkImage(img, stop); } ItkUcharImgType::Pointer tissue; if (!tissueFile.empty()) { MITK_INFO << "loading tissue image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(tissueFile)[0].GetPointer()); tissue = ItkUcharImgType::New(); mitk::CastToItkImage(img, tissue); } MITK_INFO << "loading additional images"; typedef itk::Image ItkFloatImgType; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(file)[0].GetPointer()); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); if (algorithm == "ProbTensor") { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkTensorImg = caster->GetOutput(); typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkTensorImg ); filter->Update(); mitk::Image::Pointer image = mitk::Image::New(); FilterType::OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); input_images.clear(); input_images.push_back(image); algorithm = "ProbODF"; sharpen_odfs = true; odf_cutoff = 0; } typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { mitk::TractographyForest::Pointer forest = dynamic_cast(mitk::IOUtil::Load(forestFile)[0].GetPointer()); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } if (algorithm == "ProbRF") handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); } else if (algorithm == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetPeakImage(itkImg); dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); dynamic_cast(handler)->SetPeakThreshold(cutoff); } else if (algorithm == "DetTensor") { handler = new mitk::TrackingHandlerTensor(); for (auto input_image : input_images) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_image); caster->Update(); - mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkImg = caster->GetOutput(); + mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(handler)->AddTensorImage(itkImg); } dynamic_cast(handler)->SetFaThreshold(cutoff); dynamic_cast(handler)->SetF(tend_f); dynamic_cast(handler)->SetG(tend_g); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF") { handler = new mitk::TrackingHandlerOdf(); typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetOdfImage(itkImg); dynamic_cast(handler)->SetGfaThreshold(cutoff); dynamic_cast(handler)->SetOdfThreshold(odf_cutoff); dynamic_cast(handler)->SetSharpenOdfs(sharpen_odfs); if (algorithm == "ProbODF") dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } handler->SetInterpolate(interpolate); handler->SetFlipX(flip_x); handler->SetFlipY(flip_y); handler->SetFlipZ(flip_z); MITK_INFO << "Tractography algorithm: " << algorithm; tracker->SetNumberOfSamples(num_samples); tracker->SetAngularThreshold(angular_threshold); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetSeedsPerVoxel(seeds); tracker->SetStepSize(stepsize); tracker->SetSamplingDistance(sampling_distance); tracker->SetUseStopVotes(use_stop_votes); tracker->SetOnlyForwardSamples(use_only_forward_samples); tracker->SetAposterioriCurvCheck(false); tracker->SetTissueImage(tissue); tracker->SetSeedOnlyGm(seed_gm); tracker->SetControlGmEndings(control_gm_endings); tracker->SetMaxNumTracts(max_tracts); tracker->SetTrackingHandler(handler); tracker->SetUseOutputProbabilityMap(output_prob_map); tracker->SetMinTractLength(min_tract_length); tracker->Update(); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (!output_prob_map) { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); if (ext != ".fib" && ext != ".trk" && ext != ".tck") outFile += ".fib"; mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index 197522c075..3b85ee12d7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1456 +1,1533 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkFiberProcessingView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include // ITK #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() - : QmitkAbstractView() - , m_Controls( 0 ) - , m_CircleCounter(0) - , m_PolygonCounter(0) - , m_UpsamplingFactor(1) + : QmitkAbstractView() + , m_Controls( 0 ) + , m_CircleCounter(0) + , m_PolygonCounter(0) + , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { - + RemoveObservers(); } void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent ) { - // build up qt view, unless already done - if ( !m_Controls ) - { - // create GUI widgets from the Qt Designer's .ui file - m_Controls = new Ui::QmitkFiberProcessingViewControls; - m_Controls->setupUi( parent ); - - connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); - connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); - connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); - connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); - connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); - connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); - - connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); - connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); - connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); - - connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); - connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); - connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); - - connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); - connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); - connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); - connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); - - m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); - mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); - mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); - mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); - mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); - mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); - isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); - mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); - mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); - m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); - - m_Controls->label_17->setVisible(false); - m_Controls->m_FiberExtractionFractionBox->setVisible(false); - } - - UpdateGui(); + // build up qt view, unless already done + if ( !m_Controls ) + { + // create GUI widgets from the Qt Designer's .ui file + m_Controls = new Ui::QmitkFiberProcessingViewControls; + m_Controls->setupUi( parent ); + + connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); + connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); + connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); + connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); + connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); + connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); + + connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); + connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); + connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); + + connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); + connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); + connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); + + connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); + connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); + connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); + connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); + + m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); + mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); + mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); + mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); + mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); + mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); + isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); + mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); + mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); + m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); + + m_Controls->label_17->setVisible(false); + m_Controls->m_FiberExtractionFractionBox->setVisible(false); + } + + UpdateGui(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { - if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2) - { - m_Controls->label_17->setVisible(true); - m_Controls->m_FiberExtractionFractionBox->setVisible(true); - m_Controls->m_BothEnds->setVisible(false); - } - else - { - m_Controls->label_17->setVisible(false); - m_Controls->m_FiberExtractionFractionBox->setVisible(false); - if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) - m_Controls->m_BothEnds->setVisible(true); - } + if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2) + { + m_Controls->label_17->setVisible(true); + m_Controls->m_FiberExtractionFractionBox->setVisible(true); + m_Controls->m_BothEnds->setVisible(false); + } + else + { + m_Controls->label_17->setVisible(false); + m_Controls->m_FiberExtractionFractionBox->setVisible(false); + if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) + m_Controls->m_BothEnds->setVisible(true); + } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { - switch (m_Controls->m_ModificationMethodBox->currentIndex()) - { - case 0: - { - ResampleSelectedBundlesSpline(); - break; - } - case 1: - { - ResampleSelectedBundlesLinear(); - break; - } - case 2: - { - CompressSelectedBundles(); - break; - } - case 3: - { - DoImageColorCoding(); - break; - } - case 4: - { - MirrorFibers(); - break; - } - case 5: - { - WeightFibers(); - break; - } - case 6: - { - DoCurvatureColorCoding(); - break; - } - case 7: - { - DoWeightColorCoding(); - break; - } - } + switch (m_Controls->m_ModificationMethodBox->currentIndex()) + { + case 0: + { + ResampleSelectedBundlesSpline(); + break; + } + case 1: + { + ResampleSelectedBundlesLinear(); + break; + } + case 2: + { + CompressSelectedBundles(); + break; + } + case 3: + { + DoImageColorCoding(); + break; + } + case 4: + { + MirrorFibers(); + break; + } + case 5: + { + WeightFibers(); + break; + } + case 6: + { + DoCurvatureColorCoding(); + break; + } + case 7: + { + DoWeightColorCoding(); + break; + } + } } void QmitkFiberProcessingView::WeightFibers() { - float weight = this->m_Controls->m_BundleWeightBox->value(); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->SetFiberWeights(weight); - } + float weight = this->m_Controls->m_BundleWeightBox->value(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->SetFiberWeights(weight); + } } void QmitkFiberProcessingView::Remove() { - switch (m_Controls->m_RemovalMethodBox->currentIndex()) - { - case 0: - { - RemoveDir(); - break; - } - case 1: - { - PruneBundle(); - break; - } - case 2: + switch (m_Controls->m_RemovalMethodBox->currentIndex()) + { + case 0: + { + RemoveDir(); + break; + } + case 1: + { + PruneBundle(); + break; + } + case 2: + { + ApplyCurvatureThreshold(); + break; + } + case 3: + { + RemoveWithMask(false); + break; + } + case 4: + { + RemoveWithMask(true); + break; + } + } +} + +void QmitkFiberProcessingView::Extract() +{ + switch (m_Controls->m_ExtractionMethodBox->currentIndex()) + { + case 0: + { + ExtractWithPlanarFigure(); + break; + } + case 1: + { + switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { - ApplyCurvatureThreshold(); - break; - } - case 3: { - RemoveWithMask(false); + case 0: + ExtractWithMask(true, false); break; } - case 4: { - RemoveWithMask(true); + case 1: + ExtractWithMask(true, true); break; } - } -} - -void QmitkFiberProcessingView::Extract() -{ - switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { - case 0: - { - ExtractWithPlanarFigure(); + case 2: + ExtractWithMask(false, false); break; } - case 1: { - switch (m_Controls->m_ExtractionBoxMask->currentIndex()) - { - { - case 0: - ExtractWithMask(true, false); - break; - } - { - case 1: - ExtractWithMask(true, true); - break; - } - { - case 2: - ExtractWithMask(false, false); - break; - } - { - case 3: - ExtractWithMask(false, true); - break; - } - } + case 3: + ExtractWithMask(false, true); break; } } + break; + } + } } void QmitkFiberProcessingView::PruneBundle() { - int minLength = this->m_Controls->m_PruneFibersMinBox->value(); - int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - if (!fib->RemoveShortFibers(minLength)) - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); - else if (!fib->RemoveLongFibers(maxLength)) - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); - } - RenderingManager::GetInstance()->RequestUpdateAll(); + int minLength = this->m_Controls->m_PruneFibersMinBox->value(); + int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + if (!fib->RemoveShortFibers(minLength)) + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + else if (!fib->RemoveLongFibers(maxLength)) + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { - int angle = this->m_Controls->m_CurvSpinBox->value(); - int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); - std::vector< DataNode::Pointer > nodes = m_SelectedFB; - for (auto node : nodes) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - - itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); - filter->SetInputFiberBundle(fib); - filter->SetAngularDeviation(angle); - filter->SetDistance(dist); - filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); - filter->Update(); - mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); - if (newFib->GetNumFibers()>0) - { - newFib->ColorFibersByOrientation(); - node->SetData(newFib); - } - else - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + int angle = this->m_Controls->m_CurvSpinBox->value(); + int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); + std::vector< DataNode::Pointer > nodes = m_SelectedFB; + for (auto node : nodes) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + + itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); + filter->SetInputFiberBundle(fib); + filter->SetAngularDeviation(angle); + filter->SetDistance(dist); + filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); + filter->Update(); + mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); + if (newFib->GetNumFibers()>0) + { + newFib->ColorFibersByOrientation(); + node->SetData(newFib); } - RenderingManager::GetInstance()->RequestUpdateAll(); + else + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - vnl_vector_fixed dir; - dir[0] = m_Controls->m_ExtractDirX->value(); - dir[1] = m_Controls->m_ExtractDirY->value(); - dir[2] = m_Controls->m_ExtractDirZ->value(); - fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*M_PI/180)); - } - RenderingManager::GetInstance()->RequestUpdateAll(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + vnl_vector_fixed dir; + dir[0] = m_Controls->m_ExtractDirX->value(); + dir[1] = m_Controls->m_ExtractDirY->value(); + dir[2] = m_Controls->m_ExtractDirZ->value(); + fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*M_PI/180)); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { - if (m_MaskImageNode.IsNull()) - return; + if (m_MaskImageNode.IsNull()) + return; - mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - QString name(node->GetName().c_str()); + mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - itkUCharImageType::Pointer mask = itkUCharImageType::New(); - mitk::CastToItkImage(mitkMask, mask); - mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); - if (newFib->GetNumFibers()<=0) - { - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); - continue; - } - node->SetData(newFib); + itkUCharImageType::Pointer mask = itkUCharImageType::New(); + mitk::CastToItkImage(mitkMask, mask); + mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); + if (newFib->GetNumFibers()<=0) + { + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + continue; } - RenderingManager::GetInstance()->RequestUpdateAll(); + node->SetData(newFib); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert) { - if (m_MaskImageNode.IsNull()) - return; - - mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - QString name(node->GetName().c_str()); - - itkUCharImageType::Pointer mask = itkUCharImageType::New(); - mitk::CastToItkImage(mitkMask, mask); - mitk::FiberBundle::Pointer newFib = fib->ExtractFiberSubset(mask, !onlyEnds, invert, m_Controls->m_BothEnds->isChecked(), m_Controls->m_FiberExtractionFractionBox->value()); - if (newFib->GetNumFibers()<=0) - { - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); - continue; - } - - DataNode::Pointer newNode = DataNode::New(); - newNode->SetData(newFib); + if (m_MaskImageNode.IsNull()) + return; - if (invert) - { - name += "_not"; - - if (onlyEnds) - name += "-ending-in-mask"; - else - name += "-passing-mask"; - } - else - { - if (onlyEnds) - name += "_ending-in-mask"; - else - name += "_passing-mask"; - } + mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + QString name(node->GetName().c_str()); - newNode->SetName(name.toStdString()); - GetDataStorage()->Add(newNode); - node->SetVisibility(false); + itkUCharImageType::Pointer mask = itkUCharImageType::New(); + mitk::CastToItkImage(mitkMask, mask); + mitk::FiberBundle::Pointer newFib = fib->ExtractFiberSubset(mask, !onlyEnds, invert, m_Controls->m_BothEnds->isChecked(), m_Controls->m_FiberExtractionFractionBox->value()); + if (newFib->GetNumFibers()<=0) + { + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + continue; } -} -void QmitkFiberProcessingView::GenerateRoiImage() -{ - if (m_SelectedPF.empty()) - return; + DataNode::Pointer newNode = DataNode::New(); + newNode->SetData(newFib); - mitk::BaseGeometry::Pointer geometry; - if (!m_SelectedFB.empty()) + if (invert) { - mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); - geometry = fib->GetGeometry(); + name += "_not"; + + if (onlyEnds) + name += "-ending-in-mask"; + else + name += "-passing-mask"; } - else if (m_SelectedImage) - geometry = m_SelectedImage->GetGeometry(); else - return; - - itk::Vector spacing = geometry->GetSpacing(); - spacing /= m_UpsamplingFactor; - - mitk::Point3D newOrigin = geometry->GetOrigin(); - mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); - newOrigin[0] += bounds.GetElement(0); - newOrigin[1] += bounds.GetElement(2); - newOrigin[2] += bounds.GetElement(4); - - itk::Matrix direction; - itk::ImageRegion<3> imageRegion; - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; - imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); - imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); - imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); - - m_PlanarFigureImage = itkUCharImageType::New(); - m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing - m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin - m_PlanarFigureImage->SetDirection( direction ); // Set the image direction - m_PlanarFigureImage->SetRegions( imageRegion ); - m_PlanarFigureImage->Allocate(); - m_PlanarFigureImage->FillBuffer( 0 ); - - Image::Pointer tmpImage = Image::New(); - tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); - tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); - - std::string name = m_SelectedPF.at(0)->GetName(); - WritePfToImage(m_SelectedPF.at(0), tmpImage); - for (unsigned int i=1; iGetName(); - WritePfToImage(m_SelectedPF.at(i), tmpImage); + if (onlyEnds) + name += "_ending-in-mask"; + else + name += "_passing-mask"; } - DataNode::Pointer node = DataNode::New(); - tmpImage = Image::New(); - tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); - tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); - node->SetData(tmpImage); - node->SetName(name); - this->GetDataStorage()->Add(node); + newNode->SetName(name.toStdString()); + GetDataStorage()->Add(newNode); + node->SetVisibility(false); + } +} + +void QmitkFiberProcessingView::GenerateRoiImage() +{ + if (m_SelectedPF.empty()) + return; + + mitk::BaseGeometry::Pointer geometry; + if (!m_SelectedFB.empty()) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); + geometry = fib->GetGeometry(); + } + else if (m_SelectedImage) + geometry = m_SelectedImage->GetGeometry(); + else + return; + + itk::Vector spacing = geometry->GetSpacing(); + spacing /= m_UpsamplingFactor; + + mitk::Point3D newOrigin = geometry->GetOrigin(); + mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); + newOrigin[0] += bounds.GetElement(0); + newOrigin[1] += bounds.GetElement(2); + newOrigin[2] += bounds.GetElement(4); + + itk::Matrix direction; + itk::ImageRegion<3> imageRegion; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; + imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); + imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); + imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); + + m_PlanarFigureImage = itkUCharImageType::New(); + m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing + m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin + m_PlanarFigureImage->SetDirection( direction ); // Set the image direction + m_PlanarFigureImage->SetRegions( imageRegion ); + m_PlanarFigureImage->Allocate(); + m_PlanarFigureImage->FillBuffer( 0 ); + + Image::Pointer tmpImage = Image::New(); + tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); + tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); + + std::string name = m_SelectedPF.at(0)->GetName(); + WritePfToImage(m_SelectedPF.at(0), tmpImage); + for (unsigned int i=1; iGetName(); + WritePfToImage(m_SelectedPF.at(i), tmpImage); + } + + DataNode::Pointer node = DataNode::New(); + tmpImage = Image::New(); + tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); + tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); + node->SetData(tmpImage); + node->SetName(name); + this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { - if (dynamic_cast(node->GetData())) - { - m_PlanarFigure = dynamic_cast(node->GetData()); - AccessFixedDimensionByItk_2( - image, - InternalReorientImagePlane, 3, - m_PlanarFigure->GetGeometry(), -1); - - AccessFixedDimensionByItk_2( - m_InternalImage, - InternalCalculateMaskFromPlanarFigure, - 3, 2, node->GetName() ); - } - else if (dynamic_cast(node->GetData())) - { - DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); - for (unsigned int i=0; iSize(); i++) - { - WritePfToImage(children->at(i), image); - } - } + if (dynamic_cast(node->GetData())) + { + m_PlanarFigure = dynamic_cast(node->GetData()); + AccessFixedDimensionByItk_2( + image, + InternalReorientImagePlane, 3, + m_PlanarFigure->GetGeometry(), -1); + + AccessFixedDimensionByItk_2( + m_InternalImage, + InternalCalculateMaskFromPlanarFigure, + 3, 2, node->GetName() ); + } + else if (dynamic_cast(node->GetData())) + { + DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); + for (unsigned int i=0; iSize(); i++) + { + WritePfToImage(children->at(i), image); + } + } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< float, VImageDimension > FloatImageType; - - typedef itk::ResampleImageFilter ResamplerType; - typename ResamplerType::Pointer resampler = ResamplerType::New(); - - mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); - - float upsamp = m_UpsamplingFactor; - float gausssigma = 0.5; - - // Spacing - typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); - spacing[0] = image->GetSpacing()[0] / upsamp; - spacing[1] = image->GetSpacing()[1] / upsamp; - spacing[2] = image->GetSpacing()[2]; - resampler->SetOutputSpacing( spacing ); - - // Size - typename ResamplerType::SizeType size; - size[0] = planegeo->GetExtentInMM(0) / spacing[0]; - size[1] = planegeo->GetExtentInMM(1) / spacing[1]; - size[2] = 1; - resampler->SetSize( size ); - - // Origin - typename mitk::Point3D orig = planegeo->GetOrigin(); - typename mitk::Point3D corrorig; - planegeo3D->WorldToIndex(orig,corrorig); - corrorig[0] += 0.5/upsamp; - corrorig[1] += 0.5/upsamp; - corrorig[2] += 0; - planegeo3D->IndexToWorld(corrorig,corrorig); - resampler->SetOutputOrigin(corrorig ); - - // Direction - typename ResamplerType::DirectionType direction; - typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); - for(unsigned int c=0; cSetOutputDirection( direction ); - - // Gaussian interpolation - if(gausssigma != 0) - { - double sigma[3]; - for( unsigned int d = 0; d < 3; d++ ) - sigma[d] = gausssigma * image->GetSpacing()[d]; - double alpha = 2.0; - typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; - typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); - interpolator->SetInputImage( image ); - interpolator->SetParameters( sigma, alpha ); - resampler->SetInterpolator( interpolator ); - } - else - { - typedef typename itk::LinearInterpolateImageFunction InterpolatorType; - typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); - interpolator->SetInputImage( image ); - resampler->SetInterpolator( interpolator ); - } - - resampler->SetInput( image ); - resampler->SetDefaultPixelValue(0); - resampler->Update(); - - if(additionalIndex < 0) - { - this->m_InternalImage = mitk::Image::New(); - this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); - this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); - } + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< float, VImageDimension > FloatImageType; + + typedef itk::ResampleImageFilter ResamplerType; + typename ResamplerType::Pointer resampler = ResamplerType::New(); + + mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); + + float upsamp = m_UpsamplingFactor; + float gausssigma = 0.5; + + // Spacing + typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); + spacing[0] = image->GetSpacing()[0] / upsamp; + spacing[1] = image->GetSpacing()[1] / upsamp; + spacing[2] = image->GetSpacing()[2]; + resampler->SetOutputSpacing( spacing ); + + // Size + typename ResamplerType::SizeType size; + size[0] = planegeo->GetExtentInMM(0) / spacing[0]; + size[1] = planegeo->GetExtentInMM(1) / spacing[1]; + size[2] = 1; + resampler->SetSize( size ); + + // Origin + typename mitk::Point3D orig = planegeo->GetOrigin(); + typename mitk::Point3D corrorig; + planegeo3D->WorldToIndex(orig,corrorig); + corrorig[0] += 0.5/upsamp; + corrorig[1] += 0.5/upsamp; + corrorig[2] += 0; + planegeo3D->IndexToWorld(corrorig,corrorig); + resampler->SetOutputOrigin(corrorig ); + + // Direction + typename ResamplerType::DirectionType direction; + typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); + for(unsigned int c=0; cSetOutputDirection( direction ); + + // Gaussian interpolation + if(gausssigma != 0) + { + double sigma[3]; + for( unsigned int d = 0; d < 3; d++ ) + sigma[d] = gausssigma * image->GetSpacing()[d]; + double alpha = 2.0; + typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; + typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); + interpolator->SetInputImage( image ); + interpolator->SetParameters( sigma, alpha ); + resampler->SetInterpolator( interpolator ); + } + else + { + typedef typename itk::LinearInterpolateImageFunction InterpolatorType; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + interpolator->SetInputImage( image ); + resampler->SetInterpolator( interpolator ); + } + + resampler->SetInput( image ); + resampler->SetDefaultPixelValue(0); + resampler->Update(); + + if(additionalIndex < 0) + { + this->m_InternalImage = mitk::Image::New(); + this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); + this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); + } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string ) { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::CastImageFilter< ImageType, itkUCharImageType > CastFilterType; - - // Generate mask image as new image with same header as input image and - // initialize with "1". - itkUCharImageType::Pointer newMaskImage = itkUCharImageType::New(); - newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing - newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin - newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction - newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); - newMaskImage->Allocate(); - newMaskImage->FillBuffer( 1 ); - - // Generate VTK polygon from (closed) PlanarFigure polyline - // (The polyline points are shifted by -0.5 in z-direction to make sure - // that the extrusion filter, which afterwards elevates all points by +0.5 - // in z-direction, creates a 3D object which is cut by the the plane z=0) - const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); - const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); - const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); - - vtkPolyData *polyline = vtkPolyData::New(); - polyline->Allocate( 1, 1 ); - - // Determine x- and y-dimensions depending on principal axis - int i0, i1; - switch ( axis ) - { - case 0: - i0 = 1; - i1 = 2; - break; - case 1: - i0 = 0; - i1 = 2; - break; - case 2: - default: - i0 = 0; - i1 = 1; - break; - } - - // Create VTK polydata object of polyline contour - vtkPoints *points = vtkPoints::New(); - PlanarFigure::PolyLineType::const_iterator it; - unsigned int numberOfPoints = 0; - - for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) - { - Point3D point3D; - - // Convert 2D point back to the local index coordinates of the selected image - Point2D point2D = *it; - planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); - point2D[0] -= 0.5/m_UpsamplingFactor; - point2D[1] -= 0.5/m_UpsamplingFactor; - planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); - planarFigurePlaneGeometry->Map( point2D, point3D ); - - // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil - if ( !imageGeometry3D->IsInside( point3D ) ) - { - float bounds[2] = {0,0}; - bounds[0] = - this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); - bounds[1] = - this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); - - imageGeometry3D->WorldToIndex( point3D, point3D ); - - if (point3D[i0]<0) - point3D[i0] = 0.0; - else if (point3D[i0]>bounds[0]) - point3D[i0] = bounds[0]-0.001; - - if (point3D[i1]<0) - point3D[i1] = 0.0; - else if (point3D[i1]>bounds[1]) - point3D[i1] = bounds[1]-0.001; - - points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); - numberOfPoints++; - } - else - { - imageGeometry3D->WorldToIndex( point3D, point3D ); - // Add point to polyline array - points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); - numberOfPoints++; - } + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::CastImageFilter< ImageType, itkUCharImageType > CastFilterType; + + // Generate mask image as new image with same header as input image and + // initialize with "1". + itkUCharImageType::Pointer newMaskImage = itkUCharImageType::New(); + newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing + newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin + newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction + newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); + newMaskImage->Allocate(); + newMaskImage->FillBuffer( 1 ); + + // Generate VTK polygon from (closed) PlanarFigure polyline + // (The polyline points are shifted by -0.5 in z-direction to make sure + // that the extrusion filter, which afterwards elevates all points by +0.5 + // in z-direction, creates a 3D object which is cut by the the plane z=0) + const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); + const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); + const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); + + vtkPolyData *polyline = vtkPolyData::New(); + polyline->Allocate( 1, 1 ); + + // Determine x- and y-dimensions depending on principal axis + int i0, i1; + switch ( axis ) + { + case 0: + i0 = 1; + i1 = 2; + break; + case 1: + i0 = 0; + i1 = 2; + break; + case 2: + default: + i0 = 0; + i1 = 1; + break; + } + + // Create VTK polydata object of polyline contour + vtkPoints *points = vtkPoints::New(); + PlanarFigure::PolyLineType::const_iterator it; + unsigned int numberOfPoints = 0; + + for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) + { + Point3D point3D; + + // Convert 2D point back to the local index coordinates of the selected image + Point2D point2D = *it; + planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); + point2D[0] -= 0.5/m_UpsamplingFactor; + point2D[1] -= 0.5/m_UpsamplingFactor; + planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); + planarFigurePlaneGeometry->Map( point2D, point3D ); + + // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil + if ( !imageGeometry3D->IsInside( point3D ) ) + { + float bounds[2] = {0,0}; + bounds[0] = + this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); + bounds[1] = + this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); + + imageGeometry3D->WorldToIndex( point3D, point3D ); + + if (point3D[i0]<0) + point3D[i0] = 0.0; + else if (point3D[i0]>bounds[0]) + point3D[i0] = bounds[0]-0.001; + + if (point3D[i1]<0) + point3D[i1] = 0.0; + else if (point3D[i1]>bounds[1]) + point3D[i1] = bounds[1]-0.001; + + points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); + numberOfPoints++; } - polyline->SetPoints( points ); - points->Delete(); - - vtkIdType *ptIds = new vtkIdType[numberOfPoints]; - for ( vtkIdType i = 0; i < numberOfPoints; ++i ) - ptIds[i] = i; - polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); - - // Extrude the generated contour polygon - vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); - extrudeFilter->SetInputData( polyline ); - extrudeFilter->SetScaleFactor( 1 ); - extrudeFilter->SetExtrusionTypeToNormalExtrusion(); - extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); - - // Make a stencil from the extruded polygon - vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); - polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); - - // Export from ITK to VTK (to use a VTK filter) - typedef itk::VTKImageImport< itkUCharImageType > ImageImportType; - typedef itk::VTKImageExport< itkUCharImageType > ImageExportType; - - typename ImageExportType::Pointer itkExporter = ImageExportType::New(); - itkExporter->SetInput( newMaskImage ); - - vtkImageImport *vtkImporter = vtkImageImport::New(); - this->ConnectPipelines( itkExporter, vtkImporter ); - vtkImporter->Update(); - - // Apply the generated image stencil to the input image - vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); - imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); - imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); - imageStencilFilter->ReverseStencilOff(); - imageStencilFilter->SetBackgroundValue( 0 ); - imageStencilFilter->Update(); - - // Export from VTK back to ITK - vtkImageExport *vtkExporter = vtkImageExport::New(); - vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); - vtkExporter->Update(); - - typename ImageImportType::Pointer itkImporter = ImageImportType::New(); - this->ConnectPipelines( vtkExporter, itkImporter ); - itkImporter->Update(); - - // calculate cropping bounding box - m_InternalImageMask3D = itkImporter->GetOutput(); - m_InternalImageMask3D->SetDirection(image->GetDirection()); - - itk::ImageRegionConstIterator - itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); - itk::ImageRegionIterator - itimage(image, image->GetLargestPossibleRegion()); - - itmask.GoToBegin(); - itimage.GoToBegin(); - - typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; - typename ImageType::SizeType uppersize = {{0,0,0}}; - while( !itmask.IsAtEnd() ) + else { - if(itmask.Get() == 0) - itimage.Set(0); - else - { - typename ImageType::IndexType index = itimage.GetIndex(); - typename ImageType::SizeType signedindex; - signedindex[0] = index[0]; - signedindex[1] = index[1]; - signedindex[2] = index[2]; - - lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; - lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; - lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; - - uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; - uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; - uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; - } - - ++itmask; - ++itimage; - } - - typename ImageType::IndexType index; - index[0] = lowersize[0]; - index[1] = lowersize[1]; - index[2] = lowersize[2]; - - typename ImageType::SizeType size; - size[0] = uppersize[0] - lowersize[0] + 1; - size[1] = uppersize[1] - lowersize[1] + 1; - size[2] = uppersize[2] - lowersize[2] + 1; - - itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); - - // crop internal mask - typedef itk::RegionOfInterestImageFilter< itkUCharImageType, itkUCharImageType > ROIMaskFilterType; - typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); - roi2->SetRegionOfInterest(cropRegion); - roi2->SetInput(m_InternalImageMask3D); - roi2->Update(); - m_InternalImageMask3D = roi2->GetOutput(); - - Image::Pointer tmpImage = Image::New(); - tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); - tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); - - Image::Pointer tmpImage2 = Image::New(); - tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); - const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); - - const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); - - typedef itk::ImageRegionIteratorWithIndex IteratorType; - IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); - imageIterator.GoToBegin(); - while ( !imageIterator.IsAtEnd() ) + imageGeometry3D->WorldToIndex( point3D, point3D ); + // Add point to polyline array + points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); + numberOfPoints++; + } + } + polyline->SetPoints( points ); + points->Delete(); + + vtkIdType *ptIds = new vtkIdType[numberOfPoints]; + for ( vtkIdType i = 0; i < numberOfPoints; ++i ) + ptIds[i] = i; + polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); + + // Extrude the generated contour polygon + vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); + extrudeFilter->SetInputData( polyline ); + extrudeFilter->SetScaleFactor( 1 ); + extrudeFilter->SetExtrusionTypeToNormalExtrusion(); + extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); + + // Make a stencil from the extruded polygon + vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); + polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); + + // Export from ITK to VTK (to use a VTK filter) + typedef itk::VTKImageImport< itkUCharImageType > ImageImportType; + typedef itk::VTKImageExport< itkUCharImageType > ImageExportType; + + typename ImageExportType::Pointer itkExporter = ImageExportType::New(); + itkExporter->SetInput( newMaskImage ); + + vtkImageImport *vtkImporter = vtkImageImport::New(); + this->ConnectPipelines( itkExporter, vtkImporter ); + vtkImporter->Update(); + + // Apply the generated image stencil to the input image + vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); + imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); + imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); + imageStencilFilter->ReverseStencilOff(); + imageStencilFilter->SetBackgroundValue( 0 ); + imageStencilFilter->Update(); + + // Export from VTK back to ITK + vtkImageExport *vtkExporter = vtkImageExport::New(); + vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); + vtkExporter->Update(); + + typename ImageImportType::Pointer itkImporter = ImageImportType::New(); + this->ConnectPipelines( vtkExporter, itkImporter ); + itkImporter->Update(); + + // calculate cropping bounding box + m_InternalImageMask3D = itkImporter->GetOutput(); + m_InternalImageMask3D->SetDirection(image->GetDirection()); + + itk::ImageRegionConstIterator + itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); + itk::ImageRegionIterator + itimage(image, image->GetLargestPossibleRegion()); + + itmask.GoToBegin(); + itimage.GoToBegin(); + + typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; + typename ImageType::SizeType uppersize = {{0,0,0}}; + while( !itmask.IsAtEnd() ) + { + if(itmask.Get() == 0) + itimage.Set(0); + else { - unsigned char val = imageIterator.Value(); - if (val>0) - { - itk::Index<3> index = imageIterator.GetIndex(); - Point3D point; - point[0] = index[0]; - point[1] = index[1]; - point[2] = index[2]; - - intImageGeometry3D->IndexToWorld(point, point); - pfImageGeometry3D->WorldToIndex(point, point); - - point[i0] += 0.5; - point[i1] += 0.5; - - index[0] = point[0]; - index[1] = point[1]; - index[2] = point[2]; - - if (pfImageGeometry3D->IsIndexInside(index)) - m_PlanarFigureImage->SetPixel(index, 1); - } - ++imageIterator; - } - - // Clean up VTK objects - polyline->Delete(); - extrudeFilter->Delete(); - polyDataToImageStencil->Delete(); - vtkImporter->Delete(); - imageStencilFilter->Delete(); - //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? - delete[] ptIds; + typename ImageType::IndexType index = itimage.GetIndex(); + typename ImageType::SizeType signedindex; + signedindex[0] = index[0]; + signedindex[1] = index[1]; + signedindex[2] = index[2]; + + lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; + lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; + lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; + + uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; + uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; + uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; + } + + ++itmask; + ++itimage; + } + + typename ImageType::IndexType index; + index[0] = lowersize[0]; + index[1] = lowersize[1]; + index[2] = lowersize[2]; + + typename ImageType::SizeType size; + size[0] = uppersize[0] - lowersize[0] + 1; + size[1] = uppersize[1] - lowersize[1] + 1; + size[2] = uppersize[2] - lowersize[2] + 1; + + itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); + + // crop internal mask + typedef itk::RegionOfInterestImageFilter< itkUCharImageType, itkUCharImageType > ROIMaskFilterType; + typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); + roi2->SetRegionOfInterest(cropRegion); + roi2->SetInput(m_InternalImageMask3D); + roi2->Update(); + m_InternalImageMask3D = roi2->GetOutput(); + + Image::Pointer tmpImage = Image::New(); + tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); + tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); + + Image::Pointer tmpImage2 = Image::New(); + tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); + const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); + + const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); + + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); + imageIterator.GoToBegin(); + while ( !imageIterator.IsAtEnd() ) + { + unsigned char val = imageIterator.Value(); + if (val>0) + { + itk::Index<3> index = imageIterator.GetIndex(); + Point3D point; + point[0] = index[0]; + point[1] = index[1]; + point[2] = index[2]; + + intImageGeometry3D->IndexToWorld(point, point); + pfImageGeometry3D->WorldToIndex(point, point); + + point[i0] += 0.5; + point[i1] += 0.5; + + index[0] = point[0]; + index[1] = point[1]; + index[2] = point[2]; + + if (pfImageGeometry3D->IsIndexInside(index)) + m_PlanarFigureImage->SetPixel(index, 1); + } + ++imageIterator; + } + + // Clean up VTK objects + polyline->Delete(); + extrudeFilter->Delete(); + polyDataToImageStencil->Delete(); + vtkImporter->Delete(); + imageStencilFilter->Delete(); + //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? + delete[] ptIds; } void QmitkFiberProcessingView::UpdateGui() { - m_Controls->m_FibLabel->setText("mandatory"); - m_Controls->m_PfLabel->setText("needed for extraction"); - m_Controls->m_InputData->setTitle("Please Select Input Data"); - - m_Controls->m_RemoveButton->setEnabled(false); - - m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); - m_Controls->PFCompoANDButton->setEnabled(false); - m_Controls->PFCompoORButton->setEnabled(false); - m_Controls->PFCompoNOTButton->setEnabled(false); - - m_Controls->m_GenerateRoiImage->setEnabled(false); - m_Controls->m_ExtractFibersButton->setEnabled(false); - m_Controls->m_ModifyButton->setEnabled(false); - - m_Controls->m_CopyBundle->setEnabled(false); - m_Controls->m_JoinBundles->setEnabled(false); - m_Controls->m_SubstractBundles->setEnabled(false); - - // disable alle frames - m_Controls->m_BundleWeightFrame->setVisible(false); - m_Controls->m_ExtactionFramePF->setVisible(false); - m_Controls->m_RemoveDirectionFrame->setVisible(false); - m_Controls->m_RemoveLengthFrame->setVisible(false); - m_Controls->m_RemoveCurvatureFrame->setVisible(false); - m_Controls->m_SmoothFibersFrame->setVisible(false); - m_Controls->m_CompressFibersFrame->setVisible(false); - m_Controls->m_ColorFibersFrame->setVisible(false); - m_Controls->m_MirrorFibersFrame->setVisible(false); - m_Controls->m_MaskExtractionFrame->setVisible(false); - m_Controls->m_ColorMapBox->setVisible(false); - - bool pfSelected = !m_SelectedPF.empty(); - bool fibSelected = !m_SelectedFB.empty(); - bool multipleFibsSelected = (m_SelectedFB.size()>1); - bool maskSelected = m_MaskImageNode.IsNotNull(); - bool imageSelected = m_SelectedImage.IsNotNull(); - - // toggle visibility of elements according to selected method - switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) - { - case 0: - m_Controls->m_ExtactionFramePF->setVisible(true); - break; - case 1: - m_Controls->m_MaskExtractionFrame->setVisible(true); - break; - } - - switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) - { - case 0: - m_Controls->m_RemoveDirectionFrame->setVisible(true); - if ( fibSelected ) - m_Controls->m_RemoveButton->setEnabled(true); - break; - case 1: - m_Controls->m_RemoveLengthFrame->setVisible(true); - if ( fibSelected ) - m_Controls->m_RemoveButton->setEnabled(true); - break; - case 2: - m_Controls->m_RemoveCurvatureFrame->setVisible(true); - if ( fibSelected ) - m_Controls->m_RemoveButton->setEnabled(true); - break; - case 3: - break; - case 4: - break; - } - - switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) - { - case 0: - m_Controls->m_SmoothFibersFrame->setVisible(true); - break; - case 1: - m_Controls->m_SmoothFibersFrame->setVisible(true); - break; - case 2: - m_Controls->m_CompressFibersFrame->setVisible(true); - break; - case 3: - m_Controls->m_ColorFibersFrame->setVisible(true); - m_Controls->m_ColorMapBox->setVisible(true); - break; - case 4: - m_Controls->m_MirrorFibersFrame->setVisible(true); - if (m_SelectedSurfaces.size()>0) - m_Controls->m_ModifyButton->setEnabled(true); - break; - case 5: - m_Controls->m_BundleWeightFrame->setVisible(true); - break; - case 6: - m_Controls->m_ColorFibersFrame->setVisible(true); - break; - case 7: - m_Controls->m_ColorFibersFrame->setVisible(true); - break; - } - - // are fiber bundles selected? + m_Controls->m_FibLabel->setText("mandatory"); + m_Controls->m_PfLabel->setText("needed for extraction"); + m_Controls->m_InputData->setTitle("Please Select Input Data"); + + m_Controls->m_RemoveButton->setEnabled(false); + + m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); + m_Controls->PFCompoANDButton->setEnabled(false); + m_Controls->PFCompoORButton->setEnabled(false); + m_Controls->PFCompoNOTButton->setEnabled(false); + + m_Controls->m_GenerateRoiImage->setEnabled(false); + m_Controls->m_ExtractFibersButton->setEnabled(false); + m_Controls->m_ModifyButton->setEnabled(false); + + m_Controls->m_CopyBundle->setEnabled(false); + m_Controls->m_JoinBundles->setEnabled(false); + m_Controls->m_SubstractBundles->setEnabled(false); + + // disable alle frames + m_Controls->m_BundleWeightFrame->setVisible(false); + m_Controls->m_ExtactionFramePF->setVisible(false); + m_Controls->m_RemoveDirectionFrame->setVisible(false); + m_Controls->m_RemoveLengthFrame->setVisible(false); + m_Controls->m_RemoveCurvatureFrame->setVisible(false); + m_Controls->m_SmoothFibersFrame->setVisible(false); + m_Controls->m_CompressFibersFrame->setVisible(false); + m_Controls->m_ColorFibersFrame->setVisible(false); + m_Controls->m_MirrorFibersFrame->setVisible(false); + m_Controls->m_MaskExtractionFrame->setVisible(false); + m_Controls->m_ColorMapBox->setVisible(false); + + bool pfSelected = !m_SelectedPF.empty(); + bool fibSelected = !m_SelectedFB.empty(); + bool multipleFibsSelected = (m_SelectedFB.size()>1); + bool maskSelected = m_MaskImageNode.IsNotNull(); + bool imageSelected = m_SelectedImage.IsNotNull(); + + // toggle visibility of elements according to selected method + switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) + { + case 0: + m_Controls->m_ExtactionFramePF->setVisible(true); + break; + case 1: + m_Controls->m_MaskExtractionFrame->setVisible(true); + break; + } + + switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) + { + case 0: + m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) - { - m_Controls->m_CopyBundle->setEnabled(true); - m_Controls->m_ModifyButton->setEnabled(true); - m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); - m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); - - // one bundle and one planar figure needed to extract fibers - if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) - { - m_Controls->m_InputData->setTitle("Input Data"); - m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); - m_Controls->m_ExtractFibersButton->setEnabled(true); - } - - // more than two bundles needed to join/subtract - if (multipleFibsSelected) - { - m_Controls->m_FibLabel->setText("multiple bundles selected"); - - m_Controls->m_JoinBundles->setEnabled(true); - m_Controls->m_SubstractBundles->setEnabled(true); - } - - if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) - { - m_Controls->m_InputData->setTitle("Input Data"); - m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); - m_Controls->m_ExtractFibersButton->setEnabled(true); - } - - if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) - { - m_Controls->m_InputData->setTitle("Input Data"); - m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); - m_Controls->m_RemoveButton->setEnabled(true); - } - } - - // are planar figures selected? - if (pfSelected) - { - if ( fibSelected || m_SelectedImage.IsNotNull()) - m_Controls->m_GenerateRoiImage->setEnabled(true); - - if (m_SelectedPF.size() > 1) - { - m_Controls->PFCompoANDButton->setEnabled(true); - m_Controls->PFCompoORButton->setEnabled(true); - } - else - m_Controls->PFCompoNOTButton->setEnabled(true); - } - - // is image selected - if (imageSelected || maskSelected) - { - m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); + m_Controls->m_RemoveButton->setEnabled(true); + break; + case 1: + m_Controls->m_RemoveLengthFrame->setVisible(true); + if ( fibSelected ) + m_Controls->m_RemoveButton->setEnabled(true); + break; + case 2: + m_Controls->m_RemoveCurvatureFrame->setVisible(true); + if ( fibSelected ) + m_Controls->m_RemoveButton->setEnabled(true); + break; + case 3: + break; + case 4: + break; + } + + switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) + { + case 0: + m_Controls->m_SmoothFibersFrame->setVisible(true); + break; + case 1: + m_Controls->m_SmoothFibersFrame->setVisible(true); + break; + case 2: + m_Controls->m_CompressFibersFrame->setVisible(true); + break; + case 3: + m_Controls->m_ColorFibersFrame->setVisible(true); + m_Controls->m_ColorMapBox->setVisible(true); + break; + case 4: + m_Controls->m_MirrorFibersFrame->setVisible(true); + if (m_SelectedSurfaces.size()>0) + m_Controls->m_ModifyButton->setEnabled(true); + break; + case 5: + m_Controls->m_BundleWeightFrame->setVisible(true); + break; + case 6: + m_Controls->m_ColorFibersFrame->setVisible(true); + break; + case 7: + m_Controls->m_ColorFibersFrame->setVisible(true); + break; + } + + // are fiber bundles selected? + if ( fibSelected ) + { + m_Controls->m_CopyBundle->setEnabled(true); + m_Controls->m_ModifyButton->setEnabled(true); + m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); + m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); + + // one bundle and one planar figure needed to extract fibers + if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) + { + m_Controls->m_InputData->setTitle("Input Data"); + m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); + m_Controls->m_ExtractFibersButton->setEnabled(true); + } + + // more than two bundles needed to join/subtract + if (multipleFibsSelected) + { + m_Controls->m_FibLabel->setText("multiple bundles selected"); + + m_Controls->m_JoinBundles->setEnabled(true); + m_Controls->m_SubstractBundles->setEnabled(true); + } + + if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) + { + m_Controls->m_InputData->setTitle("Input Data"); + m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); + m_Controls->m_ExtractFibersButton->setEnabled(true); + } + + if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) + { + m_Controls->m_InputData->setTitle("Input Data"); + m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); + m_Controls->m_RemoveButton->setEnabled(true); + } + } + + // are planar figures selected? + if (pfSelected) + { + if ( fibSelected || m_SelectedImage.IsNotNull()) + m_Controls->m_GenerateRoiImage->setEnabled(true); + + if (m_SelectedPF.size() > 1) + { + m_Controls->PFCompoANDButton->setEnabled(true); + m_Controls->PFCompoORButton->setEnabled(true); } + else + m_Controls->PFCompoNOTButton->setEnabled(true); + } + + // is image selected + if (imageSelected || maskSelected) + { + m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); + } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* ) { - berry::IWorkbenchPart::Pointer nullPart; - QList nodes; - OnSelectionChanged(nullPart, nodes); + berry::IWorkbenchPart::Pointer nullPart; + QList nodes; + OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* ) { + if (!m_Controls->m_InteractiveBox->isChecked()) + { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); + } +} + +void QmitkFiberProcessingView::OnStartInteraction() +{ + +} + +void QmitkFiberProcessingView::OnEndInteraction() +{ + if (m_Controls->m_InteractiveBox->isChecked()) + ExtractWithPlanarFigure(true); +} + +void QmitkFiberProcessingView::AddObservers() +{ + typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType; + for (auto node : m_SelectedPF) + { + mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); + if (figure!=nullptr) + { + // add observer for event when interaction with figure starts + SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New(); + startInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnStartInteraction); + m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand ); + + // add observer for event when interaction with figure starts + SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); + endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction); + m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); + } + } +} + +void QmitkFiberProcessingView::RemoveObservers() +{ + for (auto node : m_SelectedPF) + { + mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); + if (figure!=nullptr) + { + // figure->RemoveAllObservers(); + if (figure->HasObserver(mitk::StartInteractionPlanarFigureEvent())) + figure->RemoveObserver(m_StartInteractionObserverTag); + if (figure->HasObserver(mitk::EndInteractionPlanarFigureEvent())) + figure->RemoveObserver(m_EndInteractionObserverTag); + m_StartInteractionObserverTag = -1; + m_EndInteractionObserverTag = -1; + } + } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { - //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection - m_SelectedFB.clear(); - m_SelectedPF.clear(); - m_SelectedSurfaces.clear(); - m_SelectedImage = nullptr; - m_MaskImageNode = nullptr; - - for (auto node: nodes) - { - if ( dynamic_cast(node->GetData()) ) - m_SelectedFB.push_back(node); - else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) - m_SelectedPF.push_back(node); - else if (dynamic_cast(node->GetData())) + RemoveObservers(); + + //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection + m_SelectedFB.clear(); + m_SelectedPF.clear(); + m_SelectedSurfaces.clear(); + m_SelectedImage = nullptr; + m_MaskImageNode = nullptr; + + for (auto node: nodes) + { + if ( dynamic_cast(node->GetData()) ) + m_SelectedFB.push_back(node); + else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) + m_SelectedPF.push_back(node); + else if (dynamic_cast(node->GetData())) + { + m_SelectedImage = dynamic_cast(node->GetData()); + bool isBinary = false; + node->GetPropertyValue("binary", isBinary); + if (isBinary) + m_MaskImageNode = node; + } + else if (dynamic_cast(node->GetData())) + m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); + } + + if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) + { + int maxLayer = 0; + itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); + for (unsigned int i=0; iSize(); i++) + if (dynamic_cast(nodes->at(i)->GetData())) + { + mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); + if (sources->Size()>0) + continue; + + int layer = 0; + nodes->at(i)->GetPropertyValue("layer", layer); + if (layer>=maxLayer) { - m_SelectedImage = dynamic_cast(node->GetData()); - bool isBinary = false; - node->GetPropertyValue("binary", isBinary); - if (isBinary) - m_MaskImageNode = node; + maxLayer = layer; + m_SelectedFB.clear(); + m_SelectedFB.push_back(nodes->at(i)); } - else if (dynamic_cast(node->GetData())) - m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); - } - - if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) - { - int maxLayer = 0; - itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); - for (unsigned int i=0; iSize(); i++) - if (dynamic_cast(nodes->at(i)->GetData())) - { - mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); - if (sources->Size()>0) - continue; - - int layer = 0; - nodes->at(i)->GetPropertyValue("layer", layer); - if (layer>=maxLayer) - { - maxLayer = layer; - m_SelectedFB.clear(); - m_SelectedFB.push_back(nodes->at(i)); - } - } - } - - if (m_SelectedPF.empty()) - { - int maxLayer = 0; - itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); - for (unsigned int i=0; iSize(); i++) - if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) - { - mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); - if (sources->Size()>0) - continue; - - int layer = 0; - nodes->at(i)->GetPropertyValue("layer", layer); - if (layer>=maxLayer) - { - maxLayer = layer; - m_SelectedPF.clear(); - m_SelectedPF.push_back(nodes->at(i)); - } - } - } + } + } + + if (m_SelectedPF.empty()) + { + int maxLayer = 0; + itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); + for (unsigned int i=0; iSize(); i++) + if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) + { + mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); + if (sources->Size()>0) + continue; + + int layer = 0; + nodes->at(i)->GetPropertyValue("layer", layer); + if (layer>=maxLayer) + { + maxLayer = layer; + m_SelectedPF.clear(); + m_SelectedPF.push_back(nodes->at(i)); + } + } + } - UpdateGui(); + AddObservers(); + UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { - mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); - figure->ClosedOn(); - this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); + mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); + figure->ClosedOn(); + this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { - mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); - this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); + mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); + this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { - // initialize figure's geometry with empty geometry - mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); - figure->SetPlaneGeometry( emptygeometry ); - - //set desired data to DataNode where Planarfigure is stored - mitk::DataNode::Pointer newNode = mitk::DataNode::New(); - newNode->SetName(name.toStdString()); - newNode->SetData(figure); - newNode->SetBoolProperty("planarfigure.3drendering", true); - newNode->SetBoolProperty("planarfigure.3drendering.fill", true); - - mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); - if(figureInteractor.IsNull()) - { - figureInteractor = mitk::PlanarFigureInteractor::New(); - us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); - figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); - figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); - figureInteractor->SetDataNode(newNode); - } - - // figure drawn on the topmost layer / image - GetDataStorage()->Add(newNode ); - - for(unsigned int i = 0; i < m_SelectedPF.size(); i++) - m_SelectedPF[i]->SetSelected(false); - - newNode->SetSelected(true); - m_SelectedPF.clear(); - m_SelectedPF.push_back(newNode); - UpdateGui(); + // initialize figure's geometry with empty geometry + mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); + figure->SetPlaneGeometry( emptygeometry ); + + //set desired data to DataNode where Planarfigure is stored + mitk::DataNode::Pointer newNode = mitk::DataNode::New(); + newNode->SetName(name.toStdString()); + newNode->SetData(figure); + newNode->SetBoolProperty("planarfigure.3drendering", true); + newNode->SetBoolProperty("planarfigure.3drendering.fill", true); + + mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); + if(figureInteractor.IsNull()) + { + figureInteractor = mitk::PlanarFigureInteractor::New(); + us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); + figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); + figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); + figureInteractor->SetDataNode(newNode); + } + + // figure drawn on the topmost layer / image + GetDataStorage()->Add(newNode ); + + for(unsigned int i = 0; i < m_SelectedPF.size(); i++) + m_SelectedPF[i]->SetSelected(false); + + newNode->SetSelected(true); + m_SelectedPF.clear(); + m_SelectedPF.push_back(newNode); + UpdateGui(); } -void QmitkFiberProcessingView::ExtractWithPlanarFigure() +void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { - if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ - QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); - return; + if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ + QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); + return; + } + + std::vector fiberBundles = m_SelectedFB; + mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); + for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); + mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); + + if (interactive && m_Controls->m_InteractiveBox->isChecked()) + { + if (extFB->GetNumFibers()<=0) + continue; + + if (m_InteractiveNode.IsNull()) + { + m_InteractiveNode = mitk::DataNode::New(); + QString name("Interactive"); + m_InteractiveNode->SetName(name.toStdString()); + GetDataStorage()->Add(m_InteractiveNode); + } + fib->SetFiberColors(255, 255, 255); + fiberBundles.at(i)->SetFloatProperty("opacity", 0.01); + fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); + m_InteractiveNode->SetData(extFB); + + if (auto renderWindowPart = this->GetRenderWindowPart()) + renderWindowPart->RequestUpdate(); } - - std::vector fiberBundles = m_SelectedFB; - mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); - for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); - - mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); - if (extFB->GetNumFibers()<=0) - { - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); - continue; - } - - mitk::DataNode::Pointer node; - node = mitk::DataNode::New(); - node->SetData(extFB); - QString name(fiberBundles.at(i)->GetName().c_str()); - name += "*"; - //name += planarFigure->GetName().c_str(); - node->SetName(name.toStdString()); - fiberBundles.at(i)->SetVisibility(false); - GetDataStorage()->Add(node); - } + if (extFB->GetNumFibers()<=0) + { + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + continue; + } + + mitk::DataNode::Pointer node; + node = mitk::DataNode::New(); + node->SetData(extFB); + QString name(fiberBundles.at(i)->GetName().c_str()); + name += "*"; + node->SetName(name.toStdString()); + fiberBundles.at(i)->SetVisibility(false); + GetDataStorage()->Add(node); + } + } } void QmitkFiberProcessingView::GenerateAndComposite() { - mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); - PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); - - mitk::DataNode::Pointer newPFCNode; - newPFCNode = mitk::DataNode::New(); - newPFCNode->SetName("AND"); - newPFCNode->SetData(PFCAnd); - - AddCompositeToDatastorage(newPFCNode, m_SelectedPF); - m_SelectedPF.clear(); - m_SelectedPF.push_back(newPFCNode); - UpdateGui(); + mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); + PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); + + mitk::DataNode::Pointer newPFCNode; + newPFCNode = mitk::DataNode::New(); + newPFCNode->SetName("AND"); + newPFCNode->SetData(PFCAnd); + + AddCompositeToDatastorage(newPFCNode, m_SelectedPF); + m_SelectedPF.clear(); + m_SelectedPF.push_back(newPFCNode); + UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { - mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); - PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); - - mitk::DataNode::Pointer newPFCNode; - newPFCNode = mitk::DataNode::New(); - newPFCNode->SetName("OR"); - newPFCNode->SetData(PFCOr); - - AddCompositeToDatastorage(newPFCNode, m_SelectedPF); - m_SelectedPF.clear(); - m_SelectedPF.push_back(newPFCNode); - UpdateGui(); + mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); + PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); + + mitk::DataNode::Pointer newPFCNode; + newPFCNode = mitk::DataNode::New(); + newPFCNode->SetName("OR"); + newPFCNode->SetData(PFCOr); + + AddCompositeToDatastorage(newPFCNode, m_SelectedPF); + m_SelectedPF.clear(); + m_SelectedPF.push_back(newPFCNode); + UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { - mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); - PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); - - mitk::DataNode::Pointer newPFCNode; - newPFCNode = mitk::DataNode::New(); - newPFCNode->SetName("NOT"); - newPFCNode->SetData(PFCNot); - - AddCompositeToDatastorage(newPFCNode, m_SelectedPF); - m_SelectedPF.clear(); - m_SelectedPF.push_back(newPFCNode); - UpdateGui(); + mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); + PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); + + mitk::DataNode::Pointer newPFCNode; + newPFCNode = mitk::DataNode::New(); + newPFCNode->SetName("NOT"); + newPFCNode->SetData(PFCNot); + + AddCompositeToDatastorage(newPFCNode, m_SelectedPF); + m_SelectedPF.clear(); + m_SelectedPF.push_back(newPFCNode); + UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { - pfc->SetSelected(true); - if (parentNode.IsNotNull()) - GetDataStorage()->Add(pfc, parentNode); - else - GetDataStorage()->Add(pfc); - - for (auto child : children) - { - if (dynamic_cast(child->GetData())) - { - mitk::DataNode::Pointer newChild; - newChild = mitk::DataNode::New(); - newChild->SetData(dynamic_cast(child->GetData())); - newChild->SetName( child->GetName() ); - newChild->SetBoolProperty("planarfigure.3drendering", true); - newChild->SetBoolProperty("planarfigure.3drendering.fill", true); - - GetDataStorage()->Add(newChild, pfc); - GetDataStorage()->Remove(child); - } - else if (dynamic_cast(child->GetData())) - { - mitk::DataNode::Pointer newChild; - newChild = mitk::DataNode::New(); - newChild->SetData(dynamic_cast(child->GetData())); - newChild->SetName( child->GetName() ); - - std::vector< mitk::DataNode::Pointer > grandChildVector; - mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); - for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) - grandChildVector.push_back(*it); - - AddCompositeToDatastorage(newChild, grandChildVector, pfc); - GetDataStorage()->Remove(child); - } - } - UpdateGui(); + pfc->SetSelected(true); + if (parentNode.IsNotNull()) + GetDataStorage()->Add(pfc, parentNode); + else + GetDataStorage()->Add(pfc); + + for (auto child : children) + { + if (dynamic_cast(child->GetData())) + { + mitk::DataNode::Pointer newChild; + newChild = mitk::DataNode::New(); + newChild->SetData(dynamic_cast(child->GetData())); + newChild->SetName( child->GetName() ); + newChild->SetBoolProperty("planarfigure.3drendering", true); + newChild->SetBoolProperty("planarfigure.3drendering.fill", true); + + GetDataStorage()->Add(newChild, pfc); + GetDataStorage()->Remove(child); + } + else if (dynamic_cast(child->GetData())) + { + mitk::DataNode::Pointer newChild; + newChild = mitk::DataNode::New(); + newChild->SetData(dynamic_cast(child->GetData())); + newChild->SetName( child->GetName() ); + + std::vector< mitk::DataNode::Pointer > grandChildVector; + mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); + for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) + grandChildVector.push_back(*it); + + AddCompositeToDatastorage(newChild, grandChildVector, pfc); + GetDataStorage()->Remove(child); + } + } + UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { - if ( m_SelectedFB.empty() ){ - QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); - MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; - return; - } - - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); - - node->SetVisibility(false); - QString name(""); - name += QString(m_SelectedFB.at(0)->GetName().c_str()); - name += "_copy"; - - mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); - fbNode->SetData(newFib); - fbNode->SetName(name.toStdString()); - fbNode->SetVisibility(true); - GetDataStorage()->Add(fbNode); - } - UpdateGui(); -} - -void QmitkFiberProcessingView::JoinBundles() -{ - if ( m_SelectedFB.size()<2 ){ - QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); - MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; - return; - } - - mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); - m_SelectedFB.at(0)->SetVisibility(false); + if ( m_SelectedFB.empty() ){ + QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); + MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; + return; + } + + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); + + node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); - for (unsigned int i=1; iAddBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); - name += "+"+QString(m_SelectedFB.at(i)->GetName().c_str()); - m_SelectedFB.at(i)->SetVisibility(false); - } + name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); - fbNode->SetData(newBundle); + fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); - UpdateGui(); + } + UpdateGui(); } -void QmitkFiberProcessingView::SubstractBundles() +void QmitkFiberProcessingView::JoinBundles() { - if ( m_SelectedFB.size()<2 ){ - QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); - MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; - return; - } + if ( m_SelectedFB.size()<2 ){ + QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); + MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; + return; + } + + mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); + m_SelectedFB.at(0)->SetVisibility(false); + QString name(""); + name += QString(m_SelectedFB.at(0)->GetName().c_str()); + for (unsigned int i=1; iAddBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); + name += "+"+QString(m_SelectedFB.at(i)->GetName().c_str()); + m_SelectedFB.at(i)->SetVisibility(false); + } + + mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); + fbNode->SetData(newBundle); + fbNode->SetName(name.toStdString()); + fbNode->SetVisibility(true); + GetDataStorage()->Add(fbNode); + UpdateGui(); +} - mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); - m_SelectedFB.at(0)->SetVisibility(false); - QString name(""); - name += QString(m_SelectedFB.at(0)->GetName().c_str()); - for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); - if (newBundle.IsNull()) - break; - name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); - m_SelectedFB.at(i)->SetVisibility(false); - } +void QmitkFiberProcessingView::SubstractBundles() +{ + if ( m_SelectedFB.size()<2 ){ + QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); + MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; + return; + } + + mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); + m_SelectedFB.at(0)->SetVisibility(false); + QString name(""); + name += QString(m_SelectedFB.at(0)->GetName().c_str()); + for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) - { - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); - return; - } - - mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); - fbNode->SetData(newBundle); - fbNode->SetName(name.toStdString()); - fbNode->SetVisibility(true); - GetDataStorage()->Add(fbNode); - UpdateGui(); + break; + name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); + m_SelectedFB.at(i)->SetVisibility(false); + } + if (newBundle.IsNull()) + { + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); + return; + } + + mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); + fbNode->SetData(newBundle); + fbNode->SetName(name.toStdString()); + fbNode->SetVisibility(true); + GetDataStorage()->Add(fbNode); + UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { - double factor = this->m_Controls->m_SmoothFibersBox->value(); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ResampleSpline(factor); - } - RenderingManager::GetInstance()->RequestUpdateAll(); + double factor = this->m_Controls->m_SmoothFibersBox->value(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->ResampleSpline(factor); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { - double factor = this->m_Controls->m_SmoothFibersBox->value(); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ResampleLinear(factor); - } - RenderingManager::GetInstance()->RequestUpdateAll(); + double factor = this->m_Controls->m_SmoothFibersBox->value(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->ResampleLinear(factor); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { - double factor = this->m_Controls->m_ErrorThresholdBox->value(); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->Compress(factor); - fib->ColorFibersByOrientation(); - } - RenderingManager::GetInstance()->RequestUpdateAll(); + double factor = this->m_Controls->m_ErrorThresholdBox->value(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->Compress(factor); + fib->ColorFibersByOrientation(); + } + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { - if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) - { - QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); - return; - } - - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); - } - - if (auto renderWindowPart = this->GetRenderWindowPart()) - { - renderWindowPart->RequestUpdate(); - } + if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) + { + QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); + return; + } + + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); + } + + if (auto renderWindowPart = this->GetRenderWindowPart()) + { + renderWindowPart->RequestUpdate(); + } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); - } - - if (auto renderWindowPart = this->GetRenderWindowPart()) - { - renderWindowPart->RequestUpdate(); - } + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); + } + + if (auto renderWindowPart = this->GetRenderWindowPart()) + { + renderWindowPart->RequestUpdate(); + } } void QmitkFiberProcessingView::DoWeightColorCoding() { - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); - } - - if (auto renderWindowPart = this->GetRenderWindowPart()) - { - renderWindowPart->RequestUpdate(); - } + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); + } + + if (auto renderWindowPart = this->GetRenderWindowPart()) + { + renderWindowPart->RequestUpdate(); + } } void QmitkFiberProcessingView::MirrorFibers() { - unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); - for (auto node : m_SelectedFB) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - if (m_SelectedImage.IsNotNull()) - fib->SetReferenceGeometry(m_SelectedImage->GetGeometry()); - fib->MirrorFibers(axis); - } - - - for (auto surf : m_SelectedSurfaces) - { - vtkSmartPointer poly = surf->GetVtkPolyData(); - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - - for (int i=0; iGetNumberOfPoints(); i++) - { - double* point = poly->GetPoint(i); - point[axis] *= -1; - vtkNewPoints->InsertNextPoint(point); - } - poly->SetPoints(vtkNewPoints); - surf->CalculateBoundingBox(); - } - - if (auto renderWindowPart = this->GetRenderWindowPart()) - { - renderWindowPart->RequestUpdate(); - } + unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + if (m_SelectedImage.IsNotNull()) + fib->SetReferenceGeometry(m_SelectedImage->GetGeometry()); + fib->MirrorFibers(axis); + } + + + for (auto surf : m_SelectedSurfaces) + { + vtkSmartPointer poly = surf->GetVtkPolyData(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + + for (int i=0; iGetNumberOfPoints(); i++) + { + double* point = poly->GetPoint(i); + point[axis] *= -1; + vtkNewPoints->InsertNextPoint(point); + } + poly->SetPoints(vtkNewPoints); + surf->CalculateBoundingBox(); + } + + if (auto renderWindowPart = this->GetRenderWindowPart()) + { + renderWindowPart->RequestUpdate(); + } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h index 7c83fc44ae..b9aaaaf09e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h @@ -1,191 +1,200 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkFiberProcessingView_h #define QmitkFiberProcessingView_h #include #include "ui_QmitkFiberProcessingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /*! \brief View to process fiber bundles. Supplies methods to extract fibers from the bundle, fiber resampling, mirroring, join and subtract bundles and much more. */ class QmitkFiberProcessingView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: typedef itk::Image< unsigned char, 3 > itkUCharImageType; static const std::string VIEW_ID; QmitkFiberProcessingView(); virtual ~QmitkFiberProcessingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void OnDrawCircle(); ///< add circle interactors etc. void OnDrawPolygon(); ///< add circle interactors etc. void GenerateAndComposite(); void GenerateOrComposite(); void GenerateNotComposite(); void CopyBundles(); ///< add copies of selected bundles to data storage void JoinBundles(); ///< merge selected fiber bundles void SubstractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection. void GenerateRoiImage(); ///< generate binary image of selected planar figures. void Remove(); void Extract(); void Modify(); void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection void OnMaskExtractionChanged(); virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey = nullptr, mitk::BaseProperty *property = nullptr ); protected: void MirrorFibers(); ///< mirror bundle on the specified plane void ResampleSelectedBundlesSpline(); ///< void ResampleSelectedBundlesLinear(); ///< void DoImageColorCoding(); ///< color fibers by selected scalar image void DoWeightColorCoding(); ///< color fibers by their respective weights void DoCurvatureColorCoding(); ///< color fibers by curvature void CompressSelectedBundles(); ///< remove points below certain error threshold void WeightFibers(); void RemoveWithMask(bool removeInside); void RemoveDir(); void ApplyCurvatureThreshold(); ///< remove/split fibers with a too high curvature threshold void PruneBundle(); ///< remove too short/too long fibers void ExtractWithMask(bool onlyEnds, bool invert); - void ExtractWithPlanarFigure(); + void ExtractWithPlanarFigure(bool interactive=false); + + void OnEndInteraction(); + void OnStartInteraction(); /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkFiberProcessingViewControls* m_Controls; /** Connection from VTK to ITK */ template void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName ); template < typename TPixel, unsigned int VImageDimension > void InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ); int m_CircleCounter; ///< used for data node naming int m_PolygonCounter; ///< used for data node naming std::vector m_SelectedFB; ///< selected fiber bundle nodes std::vector m_SelectedPF; ///< selected planar figure nodes std::vector m_SelectedSurfaces; mitk::Image::Pointer m_SelectedImage; mitk::Image::Pointer m_InternalImage; mitk::PlanarFigure::Pointer m_PlanarFigure; itkUCharImageType::Pointer m_InternalImageMask3D; itkUCharImageType::Pointer m_PlanarFigureImage; float m_UpsamplingFactor; ///< upsampling factor for all image generations mitk::DataNode::Pointer m_MaskImageNode; + unsigned int m_StartInteractionObserverTag; + unsigned int m_EndInteractionObserverTag; + mitk::DataNode::Pointer m_InteractiveNode; + void AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode=nullptr); void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int ); void WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image); mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute); mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib); void NodeAdded( const mitk::DataNode* node ) override; void NodeRemoved(const mitk::DataNode* node) override; + void RemoveObservers(); + void AddObservers(); }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui index c343a7217a..d49bc8afea 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,1443 +1,1450 @@ QmitkFiberProcessingViewControls 0 0 - 386 - 540 + 385 + 564 Form 0 5 0 0 - 354 - 347 + 353 + 382 Fiber Extraction Extract a fiber subset from the selected fiber bundle using manually placed planar figures as waypoints or binary regions of interest. false 0 0 200 16777215 11 Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. Extract Qt::Vertical 20 40 QFrame::NoFrame QFrame::Raised 9 9 9 9 0 - - - - false - - - - 0 - 0 - - - - - 16777215 - 16777215 - - - - - 11 - - - - Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image. - - - Generate ROI Image - - - 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 false 60 16777215 Create OR composition with selected ROIs. OR Qt::Horizontal 40 20 false 60 16777215 Create NOT composition from selected ROI. NOT false 60 16777215 Create AND composition with selected ROIs. AND 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 30 30 Draw circular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true Qt::Horizontal 40 20 30 30 Draw polygonal ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/polygon.png:/QmitkDiffusionImaging/polygon.png 32 32 true true + + + + false + + + + 0 + 0 + + + + + 16777215 + 16777215 + + + + + 11 + + + + Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image. + + + Generate ROI Image + + + + + + + Interactive Extraction + + + 0 0 Extract using planar figures Extract using binary ROI image QFrame::NoFrame QFrame::Raised 0 0 0 0 Both ends true 1.000000000000000 0.100000000000000 Extract fibers: 0 0 Ending in mask Not ending in mask Passing mask Not passing mask Min. overlap: 0 0 - 354 - 342 + 353 + 371 Fiber Removal Remove fibers that satisfy certain criteria from the selected bundle. 0 0 Remove fibers in direction Remove fibers by length Remove fibers by curvature Remove fiber parts outside mask Remove fiber parts inside mask QFrame::NoFrame QFrame::Raised 0 0 0 0 Qt::Horizontal 40 20 Minimum fiber length in mm 0 999999999 20 Max. Length: Min. Length: Maximum fiber length in mm 0 999999999 300 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 X: Y: Z: Angle: Angular deviation threshold in degree 1 90.000000000000000 1.000000000000000 25.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 If unchecked, the fiber exceeding the threshold will be split in two instead of removed. Remove Fiber false QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Max. Angular Deviation: Qt::Horizontal 40 20 Maximum angular deviation in degree 180.000000000000000 0.100000000000000 30.000000000000000 Distance: Distance in mm 1 999.000000000000000 1.000000000000000 10.000000000000000 false 0 0 200 16777215 11 Remove Qt::Vertical 20 40 0 0 - 354 - 329 + 353 + 364 Bundle Modification Modify the selected bundle with operations such as fiber resampling, FA coloring, etc. QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Error threshold in mm: 999999999.000000000000000 0.100000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Sagittal Coronal Axial Select direction: QFrame::NoFrame QFrame::Raised 0 0 0 0 0 If checked, the image values are not only used to color the fibers but are also used as opaxity values. Values as opacity false Scalar map: The values used to color the fibers are min-max normalized. If not checked, the values should be between 0 and 1. Normalize values true 0 0 Resample fibers (spline) Resample fibers (linear) Compress fibers Color fibers by scalar map (e.g. FA) Mirror fibers Weight bundle Color fibers by curvature Color fibers by fiber weights QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0.010000000000000 999999999.000000000000000 0.100000000000000 1.000000000000000 Point distance in mm: Qt::Vertical 20 40 false 0 0 200 16777215 11 Execute QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Weight: 7 999999999.000000000000000 0.100000000000000 1.000000000000000 0 0 - 368 + 367 152 Bundle Operations Join, subtract or copy bundles. false 0 0 200 16777215 11 Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. Substract Qt::Vertical 20 40 false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Join false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Copy Please Select Input Data <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true <html><head/><body><p><span style=" color:#969696;">needed for extraction</span></p></body></html> true Input DTI Fiber Bundle: Binary seed ROI. If not specified, the whole image area is seeded. ROI: Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index 8e058b5ec6..3b2074a76a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,839 +1,842 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TissueImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestBox->SetPredicate(isTractographyForest); m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageBox->SetZeroEntryText("--"); m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->SetZeroEntryText("--"); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_TissueImageBox->SetZeroEntryText("--"); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(valueChanged(double)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(valueChanged(int)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedGmBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); StartStopTrackingGui(false); } UpdateGui(); } void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { m_TrackingTimer->stop(); if (!m_Tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_InputImageNodes.at(0)->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); GetDataStorage()->Add(node, m_InputImageNodes.at(0)); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_InputImageNodes.at(0)->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_InputImageNodes.at(0)); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if (m_ThreadIsRunning) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked())) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } void QmitkStreamlineTrackingView::OnParameterChanged() { if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedGmBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_TissueImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_10->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { int dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("mandatory"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageBox->setEnabled(false); m_Controls->mFaImageLabel->setEnabled(false); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestBox->setEnabled(false); m_Controls->m_ForestLabel->setEnabled(false); m_Controls->commandLinkButton->setEnabled(false); if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_SeedGmBox->setEnabled(true); else m_Controls->m_SeedGmBox->setEnabled(false); if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.size()+" images selected"); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->m_InputData->setTitle("Input Data"); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestBox->setEnabled(true); m_Controls->m_ForestLabel->setEnabled(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } else m_Controls->m_InputData->setTitle("Please Select Input Data"); } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_ThreadIsRunning) return; if (m_InputImages.empty()) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; typedef mitk::ImageToItk CasterType; if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } if (m_FirstTensorProbRun) { QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. ODFs are sharpened (min-max normalized and raised to the power of 4). TEND parameters are ignored."); m_FirstTensorProbRun = false; } if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); TensorImageType::Pointer itkImg = TensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { - TensorImageType::Pointer itkImg = TensorImageType::New(); - mitk::CastToItkImage(m_InputImages.at(i), itkImg); + typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(m_InputImages.at(i)); + caster->Update(); + mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(m_TrackingHandler)->AddTensorImage(itkImg); } if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestBox->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography selected."); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestBox->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); StartStopTrackingGui(false); return; } try { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(itkImg); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { StartStopTrackingGui(false); return; } } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } else if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TissueImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetTissueImage(mask); m_Tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); } m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); m_Tracker->SetAposterioriCurvCheck(false); m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h index 097653bee0..5d229f4253 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h @@ -1,143 +1,143 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkStreamlineTrackingView_h #define QmitkStreamlineTrackingView_h #include #include "ui_QmitkStreamlineTrackingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class QmitkStreamlineTrackingView; class QmitkStreamlineTrackingWorker : public QObject { Q_OBJECT public: QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view); public slots: void run(); private: QmitkStreamlineTrackingView* m_View; }; /*! \brief View for tensor based deterministic streamline fiber tracking. */ class QmitkStreamlineTrackingView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; typedef itk::StreamlineTrackingFilter TrackerType; QmitkStreamlineTrackingView(); virtual ~QmitkStreamlineTrackingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; TrackerType::Pointer m_Tracker; QmitkStreamlineTrackingWorker m_TrackingWorker; QThread m_TrackingThread; protected slots: void DoFiberTracking(); ///< start fiber tracking void UpdateGui(); void ToggleInteractive(); void DeleteTrackingHandler(); void OnParameterChanged(); void InteractiveSeedChanged(bool posChanged=false); void ForestSwitched(); void OutputStyleSwitched(); void AfterThread(); ///< update gui etc. after tracking has finished void BeforeThread(); ///< start timer etc. void TimerUpdate(); void StopTractography(); void OnSliceChanged(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkStreamlineTrackingViewControls* m_Controls; protected slots: private: void StartStopTrackingGui(bool start); std::vector< itk::Point > m_SeedPoints; mitk::DataNode::Pointer m_InteractiveNode; mitk::DataNode::Pointer m_InteractivePointSetNode; std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input image nodes - std::vector< mitk::Image::Pointer > m_InputImages; ///< input images - std::vector< mitk::Image::Pointer > m_AdditionalInputImages; + std::vector< mitk::Image::ConstPointer > m_InputImages; ///< input images + std::vector< mitk::Image::ConstPointer > m_AdditionalInputImages; bool m_FirstTensorProbRun; bool m_FirstInteractiveRun; mitk::TrackingDataHandler* m_TrackingHandler; bool m_ThreadIsRunning; QTimer* m_TrackingTimer; bool m_DeleteTrackingHandler; QmitkSliceNavigationListener m_SliceChangeListener; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED