diff --git a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h index a5f8ba883e..5ed02911f1 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h @@ -1,118 +1,118 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkDiffusionFunctionCollection_h_ #define __mitkDiffusionFunctionCollection_h_ #include #include #include #include #include #include #include namespace mitk{ class imv { public: template< class TPixelType, class TOutPixelType=TPixelType > static TOutPixelType GetImageValue(const itk::Point& itkP, itk::Image* image, bool interpolate) { itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); typename itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::Pointer interpolator = itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::New(); interpolator->SetInputImage(image); if (interpolator->IsInsideBuffer(cIdx)) { if (interpolate) return interpolator->EvaluateAtContinuousIndex(cIdx); else { itk::Index<3> idx; image->TransformPhysicalPointToIndex(itkP, idx); return image->GetPixel(idx); } } else return 0.0; } template< class TPixelType=unsigned char > - static bool IsInsideMask(const itk::Point& itkP, itk::Image* image, bool interpolate) + static bool IsInsideMask(const itk::Point& itkP, itk::Image* image, bool interpolate, float threshold=0.5) { itk::ContinuousIndex< float, 3> cIdx; image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); typename itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::Pointer interpolator = itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::New(); interpolator->SetInputImage(image); if (interpolator->IsInsideBuffer(cIdx)) { double value = 0.0; if (interpolate) value = interpolator->EvaluateAtContinuousIndex(cIdx); else { itk::Index<3> idx; image->TransformPhysicalPointToIndex(itkP, idx); value = image->GetPixel(idx); } - if (value>=0.5) + if (value>=threshold) return true; } return false; } }; class MITKDIFFUSIONCORE_EXPORT sh { public: static double factorial(int number); static void Cart2Sph(double x, double y, double z, double* cart); static double legendre0(int l); static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart); static double Yj(int m, int k, double theta, double phi); }; class MITKDIFFUSIONCORE_EXPORT gradients { private: typedef std::vector IndiciesVector; typedef std::map BValueMap; typedef itk::VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; typedef vnl_vector_fixed GradientDirectionType; public: static std::vector GetAllUniqueDirections(const BValueMap &bValueMap, GradientDirectionContainerType *refGradientsContainer ); static bool CheckForDifferingShellDirections(const BValueMap &bValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer); static vnl_matrix ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder); static vnl_matrix ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer); static mitk::gradients::GradientDirectionContainerType::Pointer CreateNormalizedUniqueGradientDirectionContainer(const BValueMap &bValueMap, const GradientDirectionContainerType * origninalGradentcontainer); }; } #endif //__mitkDiffusionFunctionCollection_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp new file mode 100644 index 0000000000..acf02c42ca --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp @@ -0,0 +1,267 @@ +/*=================================================================== + +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 __itkFiberExtractionFilter_cpp +#define __itkFiberExtractionFilter_cpp + +#include "itkFiberExtractionFilter.h" + +#define _USE_MATH_DEFINES +#include +#include +#include + +namespace itk{ + + +template< class PixelType > +FiberExtractionFilter< PixelType >::FiberExtractionFilter() + : m_DontResampleFibers(false) + , m_Mode(MODE::OVERLAP) + , m_InputType(INPUT::SCALAR_MAP) + , m_BothEnds(true) + , m_OverlapFraction(0.8) + , m_NoNegatives(false) + , m_NoPositives(false) + , m_Interpolate(false) + , m_Threshold(0.5) + , m_Labels({1}) +{ + +} + +template< class PixelType > +FiberExtractionFilter< PixelType >::~FiberExtractionFilter() +{ + +} + +template< class PixelType > +void FiberExtractionFilter< PixelType >::SetRoiImages(const std::vector &rois) +{ + for (auto roi : rois) + { + if (roi==nullptr) + { + MITK_INFO << "ROI image is NULL!"; + return; + } + } + m_RoiImages = rois; +} + +template< class PixelType > +mitk::FiberBundle::Pointer FiberExtractionFilter< PixelType >::CreateFib(std::vector< long >& ids) +{ + vtkSmartPointer weights = vtkSmartPointer::New(); + vtkSmartPointer pTmp = m_InputFiberBundle->GeneratePolyDataByIds(ids, weights); + mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); + fib->SetFiberWeights(weights); + return fib; +} + +template< class PixelType > +bool FiberExtractionFilter< PixelType >::IsPositive(const itk::Point& itkP, itk::Image* image) +{ + if( m_InputType == INPUT::SCALAR_MAP ) + return mitk::imv::IsInsideMask(itkP, image, m_Interpolate, m_Threshold); + else if( m_InputType == INPUT::LABEL_MAP ) + { + auto val = mitk::imv::GetImageValue(itkP, image, m_Interpolate); + for (auto l : m_Labels) + if (l==val) + return true; + } + else + mitkThrow() << "No valid input type selected!"; + return false; +} + +template< class PixelType > +void FiberExtractionFilter< PixelType >::ExtractOverlap(mitk::FiberBundle::Pointer fib) +{ + MITK_INFO << "Extracting fibers (min. overlap " << m_OverlapFraction << ")"; + vtkSmartPointer polydata = fib->GetFiberPolyData(); + + std::vector< std::vector< long > > positive_ids; + positive_ids.resize(m_RoiImages.size()); + + std::vector< long > negative_ids; // fibers not overlapping with ANY mask + + boost::progress_display disp(m_InputFiberBundle->GetNumFibers()); + for (int i=0; iGetNumFibers(); i++) + { + ++disp; + vtkCell* cell = polydata->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + bool positive = false; + for (unsigned int m=0; mGetPoint(j); + + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + + if ( IsPositive(itkP, roi) ) + inside++; + else + outside++; + + if ((float)inside/numPoints > m_OverlapFraction) + { + positive = true; + positive_ids[m].push_back(i); + break; + } + } + } + if (!positive) + negative_ids.push_back(i); + } + + if (!m_NoNegatives) + m_Negatives.push_back(CreateFib(negative_ids)); + if (!m_NoPositives) + for (auto ids : positive_ids) + m_Positives.push_back(CreateFib(ids)); +} + +template< class PixelType > +void FiberExtractionFilter< PixelType >::ExtractEndpoints(mitk::FiberBundle::Pointer fib) +{ + MITK_INFO << "Extracting fibers (endpoints in mask)"; + vtkSmartPointer polydata = fib->GetFiberPolyData(); + + std::vector< std::vector< long > > positive_ids; + positive_ids.resize(m_RoiImages.size()); + + std::vector< long > negative_ids; // fibers not overlapping with ANY mask + + boost::progress_display disp(m_InputFiberBundle->GetNumFibers()); + for (int i=0; iGetNumFibers(); i++) + { + ++disp; + vtkCell* cell = polydata->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + bool positive = false; + if (numPoints>1) + for (unsigned int m=0; mGetPoint(0); + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + + if ( IsPositive(itkP, roi) ) + inside++; + } + + // check second fiber point + { + double* p = points->GetPoint(numPoints-1); + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + itk::Index<3> idx; + roi->TransformPhysicalPointToIndex(itkP, idx); + if ( IsPositive(itkP, roi) ) + inside++; + } + + if (inside==2 || (inside==1 && !m_BothEnds)) + { + positive = true; + positive_ids[m].push_back(i); + } + } + if (!positive) + negative_ids.push_back(i); + } + + if (!m_NoNegatives) + m_Negatives.push_back(CreateFib(negative_ids)); + if (!m_NoPositives) + for (auto ids : positive_ids) + m_Positives.push_back(CreateFib(ids)); +} + +template< class PixelType > +void FiberExtractionFilter< PixelType >::SetLabels(const std::vector &Labels) +{ + m_Labels = Labels; +} + +template< class PixelType > +std::vector FiberExtractionFilter< PixelType >::GetNegatives() const +{ + return m_Negatives; +} + +template< class PixelType > +std::vector FiberExtractionFilter< PixelType >::GetPositives() const +{ + return m_Positives; +} + +template< class PixelType > +void FiberExtractionFilter< PixelType >::GenerateData() +{ + mitk::FiberBundle::Pointer fib = m_InputFiberBundle; + if (fib->GetNumFibers()<=0) + { + MITK_INFO << "No fibers in tractogram!"; + return; + } + + if (m_Mode==MODE::OVERLAP && !m_DontResampleFibers) + { + float minSpacing = 1; + for (auto mask : m_RoiImages) + { + for (int i=0; i<3; ++i) + if(mask->GetSpacing()[i]GetSpacing()[i]; + } + + fib = m_InputFiberBundle->GetDeepCopy(); + fib->ResampleLinear(minSpacing/5); + } + + if (m_Mode == MODE::OVERLAP) + ExtractOverlap(fib); + else if (m_Mode == MODE::ENDPOINTS) + ExtractEndpoints(fib); +} + +} + +#endif // __itkFiberExtractionFilter_cpp + + + diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h new file mode 100644 index 0000000000..c50bb64edc --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h @@ -0,0 +1,110 @@ +/*=================================================================== + +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 itkFiberExtractionFilter_h +#define itkFiberExtractionFilter_h + +// MITK +#include + +// ITK +#include + +namespace itk{ + +/** +* \brief Extract streamlines from tractograms using binary images */ + +template< class PixelType > +class FiberExtractionFilter : public ProcessObject +{ +public: + + enum MODE { + OVERLAP, + ENDPOINTS + }; + + enum INPUT { + SCALAR_MAP, ///< In this case, positive means roi image > threshold + LABEL_MAP ///< In this case, positive means roi image in labels vector + }; + + typedef FiberExtractionFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + typedef itk::Image< PixelType , 3> ItkInputImgType; + + + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + itkTypeMacro( FiberExtractionFilter, ProcessObject ) + + virtual void Update() override{ + this->GenerateData(); + } + + itkSetMacro( InputFiberBundle, mitk::FiberBundle::Pointer ) + itkSetMacro( Mode, MODE ) + itkSetMacro( InputType, INPUT ) + itkSetMacro( BothEnds, bool ) + itkSetMacro( OverlapFraction, float ) + itkSetMacro( DontResampleFibers, bool ) + itkSetMacro( NoNegatives, bool ) + itkSetMacro( NoPositives, bool ) + itkSetMacro( Interpolate, bool ) + itkSetMacro( Threshold, float ) + + void SetRoiImages(const std::vector< ItkInputImgType* > &rois); + void SetLabels(const std::vector &Labels); + + std::vector GetPositives() const; + std::vector GetNegatives() const; + +protected: + + void GenerateData() override; + + FiberExtractionFilter(); + virtual ~FiberExtractionFilter(); + + mitk::FiberBundle::Pointer CreateFib(std::vector< long >& ids); + void ExtractOverlap(mitk::FiberBundle::Pointer fib); + void ExtractEndpoints(mitk::FiberBundle::Pointer fib); + bool IsPositive(const itk::Point& itkP, itk::Image* image); + + mitk::FiberBundle::Pointer m_InputFiberBundle; + std::vector< mitk::FiberBundle::Pointer > m_Positives; + std::vector< mitk::FiberBundle::Pointer > m_Negatives; + std::vector< ItkInputImgType* > m_RoiImages; + bool m_DontResampleFibers; + MODE m_Mode; + INPUT m_InputType; + bool m_BothEnds; + float m_OverlapFraction; + bool m_NoNegatives; + bool m_NoPositives; + bool m_Interpolate; + float m_Threshold; + std::vector< unsigned short > m_Labels; +}; +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkFiberExtractionFilter.cpp" +#endif + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberExtractionTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberExtractionTest.cpp index b8b6894630..87c6c99c95 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberExtractionTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberExtractionTest.cpp @@ -1,139 +1,139 @@ /*=================================================================== 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 /**Documentation * Test if fiber transfortaiom methods work correctly */ int mitkFiberExtractionTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkFiberExtractionTest"); /// \todo Fix VTK memory leaks. Bug 18097. vtkDebugLeaks::SetExitError(0); MITK_INFO << "argc: " << argc; MITK_TEST_CONDITION_REQUIRED(argc==13,"check for input data"); omp_set_num_threads(1); try{ mitk::FiberBundle::Pointer groundTruthFibs = dynamic_cast( mitk::IOUtil::Load(argv[1]).front().GetPointer() ); mitk::FiberBundle::Pointer testFibs = dynamic_cast( mitk::IOUtil::Load(argv[2]).front().GetPointer() ); // test planar figure based extraction auto data = mitk::IOUtil::Load(argv[3])[0]; auto pf1 = mitk::DataNode::New(); pf1->SetData(data); data = mitk::IOUtil::Load(argv[4])[0]; auto pf2 = mitk::DataNode::New(); pf2->SetData(data); data = mitk::IOUtil::Load(argv[5])[0]; auto pf3 = mitk::DataNode::New(); pf3->SetData(data); mitk::StandaloneDataStorage::Pointer storage = mitk::StandaloneDataStorage::New(); mitk::PlanarFigureComposite::Pointer pfc2 = mitk::PlanarFigureComposite::New(); pfc2->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer pfcNode2 = mitk::DataNode::New(); pfcNode2->SetData(pfc2); mitk::DataStorage::SetOfObjects::Pointer set2 = mitk::DataStorage::SetOfObjects::New(); set2->push_back(pfcNode2); mitk::PlanarFigureComposite::Pointer pfc1 = mitk::PlanarFigureComposite::New(); pfc1->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer pfcNode1 = mitk::DataNode::New(); pfcNode1->SetData(pfc1); mitk::DataStorage::SetOfObjects::Pointer set1 = mitk::DataStorage::SetOfObjects::New(); set1->push_back(pfcNode1); storage->Add(pfcNode2); storage->Add(pf1, set2); storage->Add(pfcNode1, set2); storage->Add(pf2, set1); storage->Add(pf3, set1); MITK_INFO << "TEST1"; mitk::FiberBundle::Pointer extractedFibs = groundTruthFibs->ExtractFiberSubset(pfcNode2, storage); MITK_INFO << "TEST2"; MITK_TEST_CONDITION_REQUIRED(extractedFibs->Equals(testFibs),"check planar figure extraction"); MITK_INFO << "TEST3"; // test subtraction and addition mitk::FiberBundle::Pointer notExtractedFibs = groundTruthFibs->SubtractBundle(extractedFibs); MITK_INFO << argv[11]; testFibs = dynamic_cast(mitk::IOUtil::Load(argv[11]).front().GetPointer()); MITK_TEST_CONDITION_REQUIRED(notExtractedFibs->Equals(testFibs),"check bundle subtraction"); mitk::FiberBundle::Pointer joinded = extractedFibs->AddBundle(notExtractedFibs); testFibs = dynamic_cast(mitk::IOUtil::Load(argv[12]).front().GetPointer()); MITK_TEST_CONDITION_REQUIRED(joinded->Equals(testFibs),"check bundle addition"); // test binary image based extraction mitk::Image::Pointer mitkRoiImage = dynamic_cast(mitk::IOUtil::Load(argv[6]).front().GetPointer()); typedef itk::Image< unsigned char, 3 > itkUCharImageType; itkUCharImageType::Pointer itkRoiImage = itkUCharImageType::New(); mitk::CastToItkImage(mitkRoiImage, itkRoiImage); { testFibs = dynamic_cast(mitk::IOUtil::Load(argv[9]).front().GetPointer()); - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); mitk::FiberBundle::Pointer test = groundTruthFibs->GetDeepCopy(); test->ResampleLinear(0.2); extractor->SetInputFiberBundle(test); - extractor->SetMasks({itkRoiImage}); + extractor->SetRoiImages({itkRoiImage}); extractor->SetOverlapFraction(0.0); extractor->SetBothEnds(true); extractor->SetDontResampleFibers(true); - extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); + extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); extractor->Update(); mitk::FiberBundle::Pointer passing = extractor->GetPositives().at(0); MITK_TEST_CONDITION_REQUIRED(passing->Equals(testFibs),"check passing mask extraction"); } { testFibs = dynamic_cast(mitk::IOUtil::Load(argv[10]).front().GetPointer()); - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(groundTruthFibs); - extractor->SetMasks({itkRoiImage}); + extractor->SetRoiImages({itkRoiImage}); extractor->SetOverlapFraction(0.0); extractor->SetBothEnds(true); - extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); + extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); mitk::FiberBundle::Pointer ending = extractor->GetPositives().at(0); MITK_TEST_CONDITION_REQUIRED(ending->Equals(testFibs),"check ending in mask extraction"); } } catch(...) { return EXIT_FAILURE; } // always end with this! MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp index 6ba9c92202..a80edf8272 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp @@ -1,144 +1,154 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; -typedef itk::Image ItkUcharImgType; +typedef itk::Image ItkFloatImgType; -ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) +ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); - ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); - mitk::CastToItkImage(img, itkMask); - return itkMask; + ItkFloatImgType::Pointer itk_image = ItkFloatImgType::New(); + mitk::CastToItkImage(img, itk_image); + return itk_image; } /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Extraction With ROI Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Extract fibers from a tractogram using binary image ROIs"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); parser.addArgument("rois", "", mitkCommandLineParser::StringList, "ROI images:", "Images with binary ROIs", us::Any(), false); parser.addArgument("both_ends", "", mitkCommandLineParser::Bool, "Both ends:", "Fibers are extracted if both endpoints are located in the ROI.", false); parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "Extract by overlap, not by endpoints. Extract fibers that overlap to at least the provided factor (0-1) with the ROI.", -1); - parser.addArgument("invert", "", mitkCommandLineParser::Bool, "Invert:", "get streamlines not positive for any of the mask images", false); - parser.addArgument("interpolate", "", mitkCommandLineParser::Bool, "Interpolate:", "interpolate mask images", false); + parser.addArgument("invert", "", mitkCommandLineParser::Bool, "Invert:", "get streamlines not positive for any of the ROI images", false); + parser.addArgument("interpolate", "", mitkCommandLineParser::Bool, "Interpolate:", "interpolate ROI images", false); + parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image > threshold", false, 0.5); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["input"]); std::string outFib = us::any_cast(parsedArgs["out"]); mitkCommandLineParser::StringContainerType roi_files = us::any_cast(parsedArgs["rois"]); bool both_ends = false; if (parsedArgs.count("both_ends")) both_ends = us::any_cast(parsedArgs["both_ends"]); bool invert = false; if (parsedArgs.count("invert")) invert = us::any_cast(parsedArgs["invert"]); float overlap_fraction = -1; if (parsedArgs.count("overlap_fraction")) overlap_fraction = us::any_cast(parsedArgs["overlap_fraction"]); bool any_point = false; if (overlap_fraction>=0) any_point = true; bool interpolate = false; if (parsedArgs.count("interpolate")) interpolate = us::any_cast(parsedArgs["interpolate"]); + float threshold = 0.5; + if (parsedArgs.count("threshold")) + threshold = us::any_cast(parsedArgs["threshold"]); + try { // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(inFib)[0].GetPointer()); - std::vector< ItkUcharImgType::Pointer > mask_images; + std::vector< ItkFloatImgType::Pointer > roi_images; for (std::size_t i=0; i roi_images2; + for (auto roi : roi_images) + roi_images2.push_back(roi); - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(inputTractogram); - extractor->SetMasks(mask_images); + extractor->SetRoiImages(roi_images2); extractor->SetOverlapFraction(overlap_fraction); extractor->SetBothEnds(both_ends); extractor->SetInterpolate(interpolate); + extractor->SetThreshold(threshold); if (invert) extractor->SetNoPositives(true); else extractor->SetNoNegatives(true); if (!any_point) - extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); + extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); - mitk::FiberBundle::Pointer newFib; + mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(nullptr); if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib); else { - newFib->AddBundles(extractor->GetPositives()); + newFib = newFib->AddBundles(extractor->GetPositives()); mitk::IOUtil::Save(newFib, outFib); } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp index 4db8ad3a87..c6864cec9f 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp @@ -1,135 +1,135 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; -typedef itk::Image ItkUcharImgType; +typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Move Overlapping Tracts"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Move tracts that overlap with one of the reference tracts"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("reference", "r", mitkCommandLineParser::StringList, "Reference:", "reference tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output Folder:", "move input tracts that do/don't overlap here", us::Any(), false); parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "", 0.9); parser.addArgument("move_overlapping", "", mitkCommandLineParser::Bool, ":", ""); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input = us::any_cast(parsedArgs["input"]); mitkCommandLineParser::StringContainerType reference = us::any_cast(parsedArgs["reference"]); std::string out_folder = us::any_cast(parsedArgs["out"]); bool move_overlapping = false; if (parsedArgs.count("move_overlapping")) move_overlapping = us::any_cast(parsedArgs["move_overlapping"]); float overlap_threshold = 0.9; if (parsedArgs.count("overlap_fraction")) overlap_threshold = us::any_cast(parsedArgs["overlap_fraction"]); try { - itk::TractDensityImageFilter< ItkUcharImgType >::Pointer filter = itk::TractDensityImageFilter< ItkUcharImgType >::New(); + itk::TractDensityImageFilter< ItkFloatImgType >::Pointer filter = itk::TractDensityImageFilter< ItkFloatImgType >::New(); filter->SetDoFiberResampling(true); filter->SetUpsamplingFactor(0.25); filter->SetBinaryOutput(true); - std::vector< ItkUcharImgType::Pointer > masks; + std::vector< ItkFloatImgType::Pointer > masks; for (auto f : reference) { mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(f)[0].GetPointer()); filter->SetFiberBundle(fib); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect filter->Update(); masks.push_back(filter->GetOutput()); std::cout.rdbuf (old); // <-- restore } boost::progress_display disp(input.size()); for (auto f : input) { ++disp; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect bool is_overlapping = false; mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(f)[0].GetPointer()); fib->ResampleLinear(2); for (auto m : masks) { float overlap = fib->GetOverlap(m, false); if (overlap>=overlap_threshold) { is_overlapping = true; break; } } if (move_overlapping && is_overlapping) ist::CopyAFile(f, out_folder + ist::GetFilenameName(f)); else if (!move_overlapping && !is_overlapping) ist::CopyAFile(f, out_folder + ist::GetFilenameName(f)); std::cout.rdbuf (old); // <-- restore } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp index caaa87c624..2bc531941f 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/EvaluateLiFE.cpp @@ -1,248 +1,248 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; -typedef itk::Image ItkUcharImgType; -typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType; +typedef itk::Image ItkFloatImgType; +typedef std::tuple< ItkFloatImgType::Pointer, std::string > MaskType; -ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) +ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); - ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); + ItkFloatImgType::Pointer itkMask = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkMask); return itkMask; } std::vector< MaskType > get_file_list(const std::string& path) { std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(std::chrono::system_clock::now().time_since_epoch()); std::srand(ms.count()); std::vector< MaskType > mask_list; itk::Directory::Pointer dir = itk::Directory::New(); if (dir->Load(path.c_str())) { int n = dir->GetNumberOfFiles(); int num_images = 0; std::vector< int > im_indices; for (int r = 0; r < n; r++) { const char *filename = dir->GetFile(r); std::string ext = ist::GetFilenameExtension(filename); if (ext==".nii" || ext==".nii.gz" || ext==".nrrd") { ++num_images; im_indices.push_back(r); } } int c = -1; for (int r : im_indices) { c++; const char *filename = dir->GetFile(r); MITK_INFO << "Loading " << ist::GetFilenameWithoutExtension(filename); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect - MaskType m(LoadItkMaskImage(path + '/' + filename), ist::GetFilenameName(filename)); + MaskType m(LoadItkImage(path + '/' + filename), ist::GetFilenameName(filename)); mask_list.push_back(m); std::cout.rdbuf (old); // <-- restore } } return mask_list; } /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Evaluate LiFE results"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", us::Any(), false); parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", false); parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Overlap threshold used to identify true positives", 0.8); parser.addArgument("steps", "", mitkCommandLineParser::Int, "Threshold steps:", "number of weight thresholds used to calculate the ROC curve", 100); parser.addArgument("pre_filter_zeros", "", mitkCommandLineParser::Bool, "Remove zero weights:", "remove fibers with zero weights before starting the evaluation"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string fibFile = us::any_cast(parsedArgs["input"]); std::string reference_mask_folder = us::any_cast(parsedArgs["reference_mask_folder"]); std::string out_folder = us::any_cast(parsedArgs["out"]); float overlap = 0.8; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); int steps = 10; if (parsedArgs.count("steps")) steps = us::any_cast(parsedArgs["steps"]); bool pre_filter_zeros = false; if (parsedArgs.count("pre_filter_zeros")) pre_filter_zeros = us::any_cast(parsedArgs["pre_filter_zeros"]); try { std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder); if (known_tract_masks.empty()) return EXIT_FAILURE; mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); // resample fibers float minSpacing = 1; if(std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[2]) minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]; else if (std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] < std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]) minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[1]; else minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]; inputTractogram->ResampleLinear(minSpacing/5); std::vector< float > weights; for (int i=0; iGetNumFibers(); ++i) weights.push_back(inputTractogram->GetFiberWeight(i)); std::sort(weights.begin(), weights.end()); if (pre_filter_zeros) inputTractogram = inputTractogram->FilterByWeights(0.0); mitk::FiberBundle::Pointer pred_positives = inputTractogram->GetDeepCopy(); mitk::FiberBundle::Pointer pred_negatives = mitk::FiberBundle::New(nullptr); ofstream logfile; logfile.open (out_folder + "LiFE_ROC.txt"); float fpr = 1.0; float tpr = 1.0; float step = weights.back()/steps; float w = 0; if (!pre_filter_zeros) w -= step; while (pred_positives->GetNumFibers()>0 && fpr>0.001 && tpr>0.001) { w += step; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::FiberBundle::Pointer tp_tracts = mitk::FiberBundle::New(nullptr); mitk::FiberBundle::Pointer fn_tracts = mitk::FiberBundle::New(nullptr); for ( MaskType mask : known_tract_masks ) { - ItkUcharImgType::Pointer mask_image = std::get<0>(mask); + ItkFloatImgType::Pointer mask_image = std::get<0>(mask); mitk::FiberBundle::Pointer a; { - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(pred_positives); - extractor->SetMasks({mask_image}); + extractor->SetRoiImages({mask_image}); extractor->SetOverlapFraction(overlap); extractor->SetDontResampleFibers(true); - extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); + extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); extractor->Update(); a = extractor->GetPositives().at(0); } tp_tracts = tp_tracts->AddBundle(a); mitk::FiberBundle::Pointer b; { - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(pred_negatives); - extractor->SetMasks({mask_image}); + extractor->SetRoiImages({mask_image}); extractor->SetOverlapFraction(overlap); extractor->SetDontResampleFibers(true); - extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); + extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); extractor->Update(); b = extractor->GetPositives().at(0); } fn_tracts = fn_tracts->AddBundle(b); } mitk::FiberBundle::Pointer fp_tracts = pred_positives->SubtractBundle(tp_tracts); mitk::FiberBundle::Pointer tn_tracts = pred_negatives->SubtractBundle(fn_tracts); std::cout.rdbuf (old); // <-- restore float positives = tp_tracts->GetNumFibers() + fn_tracts->GetNumFibers(); float negatives = tn_tracts->GetNumFibers() + fp_tracts->GetNumFibers(); fpr = (float)fp_tracts->GetNumFibers() / negatives; tpr = (float)tp_tracts->GetNumFibers() / positives; float accuracy = 1.0; if (pred_positives->GetNumFibers()>0) accuracy = (float)tp_tracts->GetNumFibers()/pred_positives->GetNumFibers(); logfile << w << " " << fpr << " " << tpr << " " << accuracy << " \n"; MITK_INFO << "#Fibers: " << pred_positives->GetNumFibers(); MITK_INFO << "FPR/TPR: " << fpr << "/" << tpr; MITK_INFO << "Acc: " << accuracy; pred_positives = inputTractogram->FilterByWeights(w); pred_negatives = inputTractogram->FilterByWeights(w, true); } logfile.close(); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp index a4044d8b6c..0d3480e294 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractAnchorTracts.cpp @@ -1,280 +1,284 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; -typedef itk::Image ItkUcharImgType; -typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType; +typedef itk::Image ItkFloatImgType; +typedef std::tuple< ItkFloatImgType::Pointer, std::string > MaskType; void CreateFolderStructure(const std::string& path) { if (ist::PathExists(path)) ist::RemoveADirectory(path); ist::MakeDirectory(path); ist::MakeDirectory(path + "/anchor_tracts/"); ist::MakeDirectory(path + "/candidate_tracts/"); ist::MakeDirectory(path + "/implausible_tracts/"); ist::MakeDirectory(path + "/skipped_masks/"); } -ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) +ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); - ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); + ItkFloatImgType::Pointer itkMask = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkMask); return itkMask; } std::vector< MaskType > get_file_list(const std::string& path, float anchor_fraction, const std::string& skipped_path, int random_seed) { if (anchor_fraction<0) anchor_fraction = 0; else if (anchor_fraction>1.0) anchor_fraction = 1.0; std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(std::chrono::system_clock::now().time_since_epoch()); if (random_seed<0) std::srand(ms.count()); else std::srand(random_seed); MITK_INFO << "random_seed: " << random_seed; std::vector< MaskType > mask_list; itk::Directory::Pointer dir = itk::Directory::New(); int skipped = 0; if (dir->Load(path.c_str())) { int n = dir->GetNumberOfFiles(); int num_images = 0; std::vector< int > im_indices; for (int r = 0; r < n; r++) { const char *filename = dir->GetFile(r); std::string ext = ist::GetFilenameExtension(filename); if (ext==".nii" || ext==".nii.gz" || ext==".nrrd") { ++num_images; im_indices.push_back(r); } } int skipping_num = num_images * (1.0 - anchor_fraction); std::random_shuffle(im_indices.begin(), im_indices.end()); MITK_INFO << "Skipping " << skipping_num << " images"; MITK_INFO << "Number of anchors: " << num_images-skipping_num; int c = -1; for (int r : im_indices) { c++; const char *filename = dir->GetFile(r); if (c matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); matrix->Identity(); matrix->SetElement(0,0,-matrix->GetElement(0,0)); matrix->SetElement(1,1,-matrix->GetElement(1,1)); geometry->SetIndexToWorldTransformByVtkMatrix(matrix); vtkSmartPointer transformFilter = vtkSmartPointer::New(); transformFilter->SetInputData(fib->GetFiberPolyData()); transformFilter->SetTransform(geometry->GetVtkTransform()); transformFilter->Update(); mitk::FiberBundle::Pointer transformed_fib = mitk::FiberBundle::New(transformFilter->GetOutput()); return transformed_fib; } /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Extract Overlapping Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", us::Any(), false); parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", us::Any(), false); parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "GM mask:", "remove fibers not ending in the gray matter"); parser.addArgument("anchor_fraction", "", mitkCommandLineParser::Float, "Anchor fraction:", "Fraction of tracts used as anchors", 0.5); parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Overlap threshold used to identify the anchor tracts", 0.8); parser.addArgument("subsample", "", mitkCommandLineParser::Float, "Subsampling factor:", "Only use specified fraction of input fibers for the analysis", 1.0); parser.addArgument("random_seed", "", mitkCommandLineParser::Int, ":", "", -1); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string fibFile = us::any_cast(parsedArgs["input"]); std::string reference_mask_folder = us::any_cast(parsedArgs["reference_mask_folder"]); std::string out_folder = us::any_cast(parsedArgs["out"]); std::string gray_matter_mask = ""; if (parsedArgs.count("gray_matter_mask")) gray_matter_mask = us::any_cast(parsedArgs["gray_matter_mask"]); float anchor_fraction = 0.5; if (parsedArgs.count("anchor_fraction")) anchor_fraction = us::any_cast(parsedArgs["anchor_fraction"]); int random_seed = -1; if (parsedArgs.count("random_seed")) random_seed = us::any_cast(parsedArgs["random_seed"]); float overlap = 0.8; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); float subsample = 1.0; if (parsedArgs.count("subsample")) subsample = us::any_cast(parsedArgs["subsample"]); try { CreateFolderStructure(out_folder); std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder, anchor_fraction, out_folder + "/skipped_masks/", random_seed); mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); MITK_INFO << "Removing fibers not ending inside of GM"; if (gray_matter_mask.compare("")!=0) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect - ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask); + ItkFloatImgType::Pointer gm_image = LoadItkImage(gray_matter_mask); std::cout.rdbuf (old); // <-- restore - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(inputTractogram); - extractor->SetMasks({gm_image}); + extractor->SetRoiImages({gm_image}); extractor->SetBothEnds(true); - extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); + extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); mitk::FiberBundle::Pointer not_gm_fibers = extractor->GetNegatives().at(0); old = cout.rdbuf(); // <-- save std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(not_gm_fibers, out_folder + "/implausible_tracts/no_gm_endings.trk"); inputTractogram = extractor->GetPositives().at(0); std::cout.rdbuf (old); // <-- restore } std::srand(0); if (subsample<1.0) inputTractogram = inputTractogram->SubsampleFibers(subsample); mitk::FiberBundle::Pointer anchor_tractogram = mitk::FiberBundle::New(nullptr); mitk::FiberBundle::Pointer candidate_tractogram = mitk::FiberBundle::New(nullptr); if (!known_tract_masks.empty()) { MITK_INFO << "Find known tracts via overlap match"; - std::vector< ItkUcharImgType::Pointer > mask_images; + std::vector< ItkFloatImgType::Pointer > mask_images; for (auto mask : known_tract_masks) mask_images.push_back(std::get<0>(mask)); - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + std::vector< ItkFloatImgType* > roi_images2; + for (auto roi : mask_images) + roi_images2.push_back(roi); + + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(inputTractogram); - extractor->SetMasks(mask_images); + extractor->SetRoiImages(roi_images2); extractor->SetOverlapFraction(overlap); - extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); + extractor->SetMode(itk::FiberExtractionFilter::MODE::OVERLAP); extractor->Update(); std::vector< mitk::FiberBundle::Pointer > positives = extractor->GetPositives(); candidate_tractogram = extractor->GetNegatives().at(0); for ( unsigned int i=0; i(known_tract_masks.at(i)); mitk::IOUtil::Save(positives.at(i), out_folder + "/anchor_tracts/" + mask_name + ".trk"); std::cout.rdbuf (old); // <-- restore } anchor_tractogram = anchor_tractogram->AddBundles(positives); } mitk::IOUtil::Save(anchor_tractogram, out_folder + "/anchor_tracts/anchor_tractogram.trk"); mitk::IOUtil::Save(candidate_tractogram, out_folder + "/candidate_tracts/candidate_tractogram.trk"); mitk::IOUtil::Save(inputTractogram, out_folder + "/filtered_tractogram.trk"); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp index dfa2276905..fd07b7bac2 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp @@ -1,240 +1,240 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include typedef itksys::SystemTools ist; -typedef itk::Image ItkUcharImgType; +typedef itk::Image ItkFloatImgType; mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } -ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) +ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); - ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); + ItkFloatImgType::Pointer itkMask = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkMask); return itkMask; } /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Extract Similar Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("ref_tracts", "", mitkCommandLineParser::StringList, "Ref. Tracts:", "reference tracts (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("ref_masks", "", mitkCommandLineParser::StringList, "Ref. Masks:", "reference bundle masks", us::Any()); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("distance", "", mitkCommandLineParser::Int, "Distance:", "", 10); parser.addArgument("metric", "", mitkCommandLineParser::String, "Metric:", ""); parser.addArgument("subsample", "", mitkCommandLineParser::Float, "Subsampling factor:", "Only use specified fraction of input fibers", 1.0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string in_fib = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType ref_bundle_files = us::any_cast(parsedArgs["ref_tracts"]); mitkCommandLineParser::StringContainerType ref_mask_files; if (parsedArgs.count("ref_masks")) ref_mask_files = us::any_cast(parsedArgs["ref_masks"]); if (ref_mask_files.size()>0 && ref_mask_files.size()!=ref_bundle_files.size()) { MITK_INFO << "If reference masks are used, there has to be one mask per reference tract."; return EXIT_FAILURE; } int distance = 10; if (parsedArgs.count("distance")) distance = us::any_cast(parsedArgs["distance"]); std::string metric = "EU_MEAN"; if (parsedArgs.count("metric")) metric = us::any_cast(parsedArgs["metric"]); float subsample = 1.0; if (parsedArgs.count("subsample")) subsample = us::any_cast(parsedArgs["subsample"]); try { mitk::FiberBundle::Pointer fib = LoadFib(in_fib); std::srand(0); if (subsample<1.0) fib = fib->SubsampleFibers(subsample); mitk::FiberBundle::Pointer resampled_fib = fib->GetDeepCopy(); resampled_fib->ResampleToNumPoints(12); std::vector< mitk::FiberBundle::Pointer > ref_fibs; - std::vector< ItkUcharImgType::Pointer > ref_masks; + std::vector< ItkFloatImgType::Pointer > ref_masks; for (std::size_t i=0; i distances; distances.push_back(distance); mitk::FiberBundle::Pointer anchor_tractogram = mitk::FiberBundle::New(nullptr); unsigned int c = 0; for (auto ref_fib : ref_fibs) { MITK_INFO << "Extracting " << ist::GetFilenameName(ref_bundle_files.at(c)); // std::streambuf *old = cout.rdbuf(); // <-- save // std::stringstream ss; // std::cout.rdbuf (ss.rdbuf()); // <-- redirect try { itk::TractClusteringFilter::Pointer segmenter = itk::TractClusteringFilter::New(); // calculate centroids from reference bundle { MITK_INFO << "TEST 1"; itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances({10,20,30}); clusterer->SetTractogram(ref_fib); clusterer->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); clusterer->SetMergeDuplicateThreshold(0.0); clusterer->Update(); MITK_INFO << "TEST 2"; std::vector tracts = clusterer->GetOutCentroids(); ref_fib = mitk::FiberBundle::New(nullptr); ref_fib = ref_fib->AddBundles(tracts); MITK_INFO << "TEST 3"; mitk::IOUtil::Save(ref_fib, out_root + "centroids_" + ist::GetFilenameName(ref_bundle_files.at(c))); segmenter->SetInCentroids(ref_fib); MITK_INFO << "TEST 4"; } // segment tract segmenter->SetFilterMask(ref_masks.at(c)); segmenter->SetOverlapThreshold(0.8); segmenter->SetDistances(distances); segmenter->SetTractogram(resampled_fib); segmenter->SetMergeDuplicateThreshold(0.0); segmenter->SetDoResampling(false); if (metric=="EU_MEAN") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); else if (metric=="EU_STD") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); else if (metric=="EU_MAX") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMax()}); segmenter->Update(); std::vector< std::vector< long > > clusters = segmenter->GetOutFiberIndices(); if (clusters.size()>0) { vtkSmartPointer weights = vtkSmartPointer::New(); mitk::FiberBundle::Pointer result = mitk::FiberBundle::New(nullptr); std::vector< mitk::FiberBundle::Pointer > result_fibs; for (unsigned int cluster_index=0; cluster_indexGeneratePolyDataByIds(clusters.at(cluster_index), weights))); result = result->AddBundles(result_fibs); anchor_tractogram = anchor_tractogram->AddBundle(result); mitk::IOUtil::Save(result, out_root + "anchor_" + ist::GetFilenameName(ref_bundle_files.at(c))); fib = mitk::FiberBundle::New(fib->GeneratePolyDataByIds(clusters.back(), weights)); resampled_fib = mitk::FiberBundle::New(resampled_fib->GeneratePolyDataByIds(clusters.back(), weights)); } } catch(itk::ExceptionObject& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.GetDescription(); } catch(std::exception& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.what(); } // std::cout.rdbuf (old); // <-- restore if (fib->GetNumFibers()==0) break; ++c; } MITK_INFO << "Streamlines in anchor tractogram: " << anchor_tractogram->GetNumFibers(); mitk::IOUtil::Save(anchor_tractogram, out_root + "anchor_tractogram.trk"); MITK_INFO << "Streamlines remaining in candidate tractogram: " << fib->GetNumFibers(); mitk::IOUtil::Save(fib, out_root + "candidate_tractogram.trk"); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp index a9488c0b72..f95146e500 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp @@ -1,251 +1,251 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; -typedef itk::Image ItkUcharImgType; +typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUIntImgType; std::vector< std::string > get_file_list(const std::string& path, std::vector< std::string > extensions={".fib", ".trk"}) { std::vector< std::string > file_list; itk::Directory::Pointer dir = itk::Directory::New(); if (dir->Load(path.c_str())) { int n = dir->GetNumberOfFiles(); for (int r = 0; r < n; r++) { const char *filename = dir->GetFile(r); std::string ext = ist::GetFilenameExtension(filename); for (auto e : extensions) { if (ext==e) { file_list.push_back(path + '/' + filename); break; } } } } return file_list; } /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Merge Overlapping Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("in", "i", mitkCommandLineParser::InputFile, "Input Folder:", "input folder", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output Folder:", "output folder", us::Any(), false); parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Tracts with overlap larger than this threshold are merged", false, 0.8); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string input_folder = us::any_cast(parsedArgs["in"]); std::string out_folder = us::any_cast(parsedArgs["out"]); float overlap = 0.8; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); try { if (ist::PathExists(out_folder)) ist::RemoveADirectory(out_folder); ist::MakeDirectory(out_folder); std::vector< std::string > fib_files = get_file_list(input_folder, {".fib", ".trk", ".tck"}); if (fib_files.empty()) return EXIT_FAILURE; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect std::vector< mitk::FiberBundle::Pointer > fibs; for (std::string f : fib_files) { mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(f)[0].GetPointer()); fibs.push_back(fib); } mitk::FiberBundle::Pointer combined = mitk::FiberBundle::New(); combined = combined->AddBundles(fibs); - itk::TractsToFiberEndingsImageFilter< ItkUcharImgType >::Pointer endings = itk::TractsToFiberEndingsImageFilter< ItkUcharImgType >::New(); + itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::Pointer endings = itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::New(); endings->SetFiberBundle(combined); endings->SetUpsamplingFactor(0.25); endings->Update(); - ItkUcharImgType::Pointer ref_image = endings->GetOutput(); + ItkFloatImgType::Pointer ref_image = endings->GetOutput(); std::cout.rdbuf (old); // <-- restore for (int its = 0; its<3; its++) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect - std::vector< ItkUcharImgType::Pointer > mask_images; + std::vector< ItkFloatImgType::Pointer > mask_images; for (auto fib : fibs) { - itk::TractDensityImageFilter< ItkUcharImgType >::Pointer masks = itk::TractDensityImageFilter< ItkUcharImgType >::New(); + itk::TractDensityImageFilter< ItkFloatImgType >::Pointer masks = itk::TractDensityImageFilter< ItkFloatImgType >::New(); masks->SetInputImage(ref_image); masks->SetBinaryOutput(true); masks->SetFiberBundle(fib); masks->SetUseImageGeometry(true); masks->Update(); mask_images.push_back(masks->GetOutput()); } int r=0; vnl_matrix< int > mat; mat.set_size(mask_images.size(), mask_images.size()); mat.fill(0); for (auto m1 : mask_images) { float max_overlap = overlap; int c = 0; for (auto m2 : mask_images) { if (c<=r) { ++c; continue; } - itk::ImageRegionConstIterator it1(m1, m1->GetLargestPossibleRegion()); - itk::ImageRegionConstIterator it2(m2, m2->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator it1(m1, m1->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator it2(m2, m2->GetLargestPossibleRegion()); unsigned int c1 = 0; unsigned int c2 = 0; unsigned int intersect = 0; while( !it1.IsAtEnd() ) { if( it1.Get()>0 && it2.Get()>0) ++intersect; if(it1.Get()>0) ++c1; if(it2.Get()>0) ++c2; ++it1; ++it2; } if ( (float)intersect/c1>max_overlap ) { max_overlap = (float)intersect/c1; mat.put(r,c, 1); } if ( (float)intersect/c2>max_overlap ) { max_overlap = (float)intersect/c2; mat.put(r,c, 1); } ++c; } ++r; } std::vector< mitk::FiberBundle::Pointer > out_fibs; std::vector< bool > used; for (unsigned int i=0; i0) { fib = fib->AddBundle(fibs.at(c)); MITK_INFO << c; used[c] = true; } } out_fibs.push_back(fib); } std::cout.rdbuf (old); // <-- restore MITK_INFO << fibs.size() << " --> " << out_fibs.size(); if (fibs.size()==out_fibs.size()) break; fibs = out_fibs; } int c = 0; for (auto fib : fibs) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + "/bundle_" + boost::lexical_cast(c) + ".trk"); std::cout.rdbuf (old); // <-- restore ++c; } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } 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 6bc0c32908..25ba19ad35 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,1598 +1,1596 @@ /*=================================================================== 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 #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) { } // 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 isOdf = mitk::NodePredicateDataType::New("OdfImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf); 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) + if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2 || m_Controls->m_ExtractionBoxMask->currentIndex() == 3) { 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; } } } 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); } } void QmitkFiberProcessingView::Remove() { 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; } case 5: { ApplyWeightThreshold(); break; } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); 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; } } 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(); } void QmitkFiberProcessingView::ApplyWeightThreshold() { float thr = this->m_Controls->m_WeightThresholdBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->FilterByWeights(thr); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByFiberWeights(false, true); node->SetData(newFib); } else 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."); } 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(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { - if (m_MaskImageNode.IsNull()) + if (m_RoiImageNode.IsNull()) return; - mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); + mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - itkUCharImageType::Pointer mask = itkUCharImageType::New(); + 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); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert) { - if (m_MaskImageNode.IsNull()) + if (m_RoiImageNode.IsNull()) return; - mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); + mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->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(); + ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(mitkMask, mask); -// mitk::FiberBundle::Pointer newFib = fib->ExtractFiberSubset(mask, !onlyEnds, invert, m_Controls->m_BothEnds->isChecked(), m_Controls->m_FiberExtractionFractionBox->value()); - itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(fib); - extractor->SetMasks({mask}); + extractor->SetRoiImages({mask}); extractor->SetOverlapFraction(m_Controls->m_FiberExtractionFractionBox->value()); extractor->SetBothEnds(m_Controls->m_BothEnds->isChecked()); + extractor->SetInterpolate(m_Controls->m_InterpolateRoiBox->isChecked()); if (invert) extractor->SetNoPositives(true); else extractor->SetNoNegatives(true); if (onlyEnds) - extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); + extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); mitk::FiberBundle::Pointer newFib; if (invert) newFib = extractor->GetNegatives().at(0); else newFib = extractor->GetPositives().at(0); if (newFib.IsNull() || 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 (invert) { name += "_not"; if (onlyEnds) name += "-ending-in-mask"; else name += "-passing-mask"; } else { if (onlyEnds) name += "_ending-in-mask"; else name += "_passing-mask"; } 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 = 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); } } } 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() ); } } 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; + 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(); + 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++; } } 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; + 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 + 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 { 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; + 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; + 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_RemoveByWeightFrame->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 maskSelected = m_RoiImageNode.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; case 5: m_Controls->m_RemoveByWeightFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); 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_PfLabel->setText(QString(m_RoiImageNode->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_PfLabel->setText(QString(m_RoiImageNode->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* node ) { for (auto fnode: m_SelectedFB) if (node == fnode) { m_SelectedFB.clear(); break; } 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::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) { figure->RemoveAllObservers(); // 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(); } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { RemoveObservers(); //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection std::vector lastSelectedFB = m_SelectedFB; m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; - m_MaskImageNode = nullptr; + m_RoiImageNode = 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; + if (m_SelectedImage->GetDimension()==3) + m_RoiImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked()) m_SelectedFB = lastSelectedFB; // if no fibers or surfaces are selected, select topmost 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 no plar figure is selected, select topmost 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)); } } } AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { 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)); } 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 ); RemoveObservers(); 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); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } try { 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 (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } float op = 5.0/sqrt(fib->GetNumFibers()); float currentOp = 0; fiberBundles.at(i)->GetFloatProperty("opacity", currentOp); if (currentOp!=op) { fib->SetFiberColors(255, 255, 255); fiberBundles.at(i)->SetFloatProperty("opacity", op); fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); } m_InteractiveNode->SetData(extFB); } else { 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); } } } catch(const std::out_of_range& ) { QMessageBox::warning( nullptr, "Fiber extraction failed", "Did you only create the planar figure, using the circle or polygon button, but forgot to actually place it in the image afterwards? \nAfter creating a planar figure, simply left-click at the desired position in the image or on the tractogram to place it."); } } 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); RemoveObservers(); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); 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); RemoveObservers(); 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); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); 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(); } 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; } m_SelectedFB.at(0)->SetVisibility(false); mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int i=1; iSetVisibility(false); tractograms.push_back(dynamic_cast(m_SelectedFB.at(i)->GetData())); } newBundle = newBundle->AddBundles(tractograms); mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName("Joined_Tractograms"); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } 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()) 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(); } 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(); } 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(); } 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(); } } 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(); } } 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(); } } 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(); } } 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 4651a3f7ed..0e56b555cf 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,200 +1,201 @@ /*=================================================================== 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; + typedef itk::Image< unsigned char, 3 > ItkUCharImageType; + typedef itk::Image< float, 3 > ItkFloatImageType; 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 ApplyWeightThreshold(); 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(bool interactive=false); void OnEndInteraction(); /// \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; + ItkUCharImageType::Pointer m_InternalImageMask3D; + ItkUCharImageType::Pointer m_PlanarFigureImage; float m_UpsamplingFactor; ///< upsampling factor for all image generations - mitk::DataNode::Pointer m_MaskImageNode; + mitk::DataNode::Pointer m_RoiImageNode; 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 627a3e8a96..3a3acd8ae6 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,1505 +1,1515 @@ QmitkFiberProcessingViewControls 0 0 385 684 Form 0 5 0 0 - 367 - 382 + 353 + 411 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 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 + Extract using ROI image QFrame::NoFrame QFrame::Raised 0 0 0 0 - + Both ends true + + + + Extract fibers: + + + 1.000000000000000 0.100000000000000 - - + + - Extract fibers: + Min. overlap: 0 0 - Ending in mask + Ending in ROI - Not ending in mask + Not ending in ROI - Passing mask + Passing ROI - Not passing mask + Not passing ROI - - + + - Min. overlap: + Interpolate ROI + + + false 0 0 367 408 Fiber Removal Remove fibers that satisfy certain criteria from the selected bundle. 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 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 Qt::Horizontal 40 20 Minimum fiber length in mm 0 999999999 20 Max. Length: Min. Length: Maximum fiber length in mm 0 999999999 300 false 0 0 200 16777215 11 Remove Qt::Vertical 20 40 0 0 Remove fibers in direction Remove fibers by length Remove fibers by curvature Remove fiber parts outside mask Remove fiber parts inside mask Remove fibers by weight QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Weight threshold: Only fibers with weight larger than this threshold are kept. 5 99999.000000000000000 0.100000000000000 0 0 367 408 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 367 165 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