diff --git a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h index 4ded475b7c..22261eac41 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h @@ -1,126 +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. ===================================================================*/ #ifndef __mitkDiffusionFunctionCollection_h_ #define __mitkDiffusionFunctionCollection_h_ #include #include #include #include #include #include #include #include namespace mitk{ class MITKDIFFUSIONCORE_EXPORT imv { public: - static std::vector< std::pair< itk::Index<3>, double > > IntersectImage(itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef); + static std::vector< std::pair< itk::Index<3>, double > > IntersectImage(const itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef); + + static itk::Point GetItkPoint(double point[3]) + { + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + return itkPoint; + } template< class TPixelType, class TOutPixelType=TPixelType > static TOutPixelType GetImageValue(const itk::Point& itkP, bool interpolate, typename itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::Pointer interpolator) { if (interpolator==nullptr) return 0.0; itk::ContinuousIndex< float, 3> cIdx; interpolator->ConvertPointToContinuousIndex(itkP, cIdx); if (interpolator->IsInsideBuffer(cIdx)) { if (interpolate) return interpolator->EvaluateAtContinuousIndex(cIdx); else { itk::Index<3> idx; interpolator->ConvertContinuousIndexToNearestIndex(cIdx, idx); return interpolator->EvaluateAtIndex(idx); } } else return 0.0; } template< class TPixelType=unsigned char > static bool IsInsideMask(const itk::Point& itkP, bool interpolate, typename itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::Pointer interpolator, float threshold=0.5) { if (interpolator==nullptr) return false; itk::ContinuousIndex< float, 3> cIdx; interpolator->ConvertPointToContinuousIndex(itkP, cIdx); if (interpolator->IsInsideBuffer(cIdx)) { double value = 0.0; if (interpolate) value = interpolator->EvaluateAtContinuousIndex(cIdx); else { itk::Index<3> idx; interpolator->ConvertContinuousIndexToNearestIndex(cIdx, idx); value = interpolator->EvaluateAtIndex(idx); } 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* spherical); 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, float theta, float phi, bool mrtrix=true); static vnl_matrix CalcShBasisForDirections(int sh_order, vnl_matrix U, bool mrtrix=true); }; class MITKDIFFUSIONCORE_EXPORT gradients { private: typedef std::vector IndiciesVector; typedef mitk::BValueMapProperty::BValueMap BValueMap; typedef DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; public: static GradientDirectionContainerType::Pointer ReadBvalsBvecs(std::string bvals_file, std::string bvecs_file, double& reference_bval); static void WriteBvalsBvecs(std::string bvals_file, std::string bvecs_file, GradientDirectionContainerType::Pointer gradients, double reference_bval); 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 GradientDirectionContainerType::Pointer CreateNormalizedUniqueGradientDirectionContainer(const BValueMap &bValueMap, const GradientDirectionContainerType * origninalGradentcontainer); }; } #endif //__mitkDiffusionFunctionCollection_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp index cf2f32da62..db526ca271 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp @@ -1,508 +1,508 @@ /*=================================================================== 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 "mitkDiffusionFunctionCollection.h" #include "mitkNumericTypes.h" #include #include #include #include #include "itkVectorContainer.h" #include "vnl/vnl_vector.h" #include #include #include #include // Intersect a finite line (with end points p0 and p1) with all of the // cells of a vtkImageData -std::vector< std::pair< itk::Index<3>, double > > mitk::imv::IntersectImage(itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef) +std::vector< std::pair< itk::Index<3>, double > > mitk::imv::IntersectImage(const itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef) { std::vector< std::pair< itk::Index<3>, double > > out; if (si == ei) { double d[3]; for (int i=0; i<3; ++i) d[i] = (sf[i]-ef[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); out.push_back( std::pair< itk::Index<3>, double >(si, len) ); return out; } double bounds[6]; double entrancePoint[3]; double exitPoint[3]; double startPoint[3]; double endPoint[3]; double t0, t1; for (int i=0; i<3; ++i) { startPoint[i] = sf[i]; endPoint[i] = ef[i]; if (si[i]>ei[i]) { int t = si[i]; si[i] = ei[i]; ei[i] = t; } } for (int x = si[0]; x<=ei[0]; ++x) for (int y = si[1]; y<=ei[1]; ++y) for (int z = si[2]; z<=ei[2]; ++z) { bounds[0] = (double)x - 0.5; bounds[1] = (double)x + 0.5; bounds[2] = (double)y - 0.5; bounds[3] = (double)y + 0.5; bounds[4] = (double)z - 0.5; bounds[5] = (double)z + 0.5; int entryPlane; int exitPlane; int hit = vtkBox::IntersectWithLine(bounds, startPoint, endPoint, t0, t1, entrancePoint, exitPoint, entryPlane, exitPlane); if (hit) { if (entryPlane>=0 && exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i] - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (entryPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (ef[i] - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i]-sf[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } } } return out; } //------------------------- SH-function ------------------------------------ double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } void mitk::sh::Cart2Sph(double x, double y, double z, double *spherical) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); if( rad < mitk::eps ) { th = itk::Math::pi/2; phi = itk::Math::pi/2; } else { th = acos(z/rad); phi = atan2(y, x); } spherical[0] = phi; spherical[1] = th; spherical[2] = rad; } double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i(l,abs(m),-cos(theta)); double mag = sqrt((double)(2*l+1)/(4.0*itk::Math::pi)*::boost::math::factorial(l-abs(m))/::boost::math::factorial(l+abs(m)))*plm; if (m>0) return mag*cos(m*phi); else if (m==0) return mag; else return mag*sin(-m*phi); } return 0; } vnl_matrix mitk::sh::CalcShBasisForDirections(int sh_order, vnl_matrix U, bool mrtrix) { vnl_matrix sh_basis = vnl_matrix(U.cols(), (sh_order*sh_order + sh_order + 2)/2 + sh_order ); for(unsigned int i=0; i bvec_entries; if (!itksys::SystemTools::FileExists(bvecs_file)) mitkThrow() << "bvecs file not existing: " << bvecs_file; else { std::string line; std::ifstream myfile (bvecs_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; boost::split(strs,line,boost::is_any_of("\t \n")); for (auto token : strs) { if (!token.empty()) { try { bvec_entries.push_back(boost::lexical_cast(token)); } catch(...) { mitkThrow() << "Encountered invalid bvecs file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvecs file could not be opened: " << bvals_file; } } reference_bval = -1; std::vector bval_entries; if (!itksys::SystemTools::FileExists(bvals_file)) mitkThrow() << "bvals file not existing: " << bvals_file; else { std::string line; std::ifstream myfile (bvals_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; boost::split(strs,line,boost::is_any_of("\t \n")); for (auto token : strs) { if (!token.empty()) { try { bval_entries.push_back(boost::lexical_cast(token)); if (bval_entries.back()>reference_bval) reference_bval = bval_entries.back(); } catch(...) { mitkThrow() << "Encountered invalid bvals file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvals file could not be opened: " << bvals_file; } } for(unsigned int i=0; i 0) { vec.normalize(); vec[0] = sqrt(factor)*vec[0]; vec[1] = sqrt(factor)*vec[1]; vec[2] = sqrt(factor)*vec[2]; } directioncontainer->InsertElement(i,vec); } return directioncontainer; } void mitk::gradients::WriteBvalsBvecs(std::string bvals_file, std::string bvecs_file, GradientDirectionContainerType::Pointer gradients, double reference_bval) { std::ofstream myfile; myfile.open (bvals_file.c_str()); for(unsigned int i=0; iSize(); i++) { double twonorm = gradients->ElementAt(i).two_norm(); myfile << std::round(reference_bval*twonorm*twonorm) << " "; } myfile.close(); std::ofstream myfile2; myfile2.open (bvecs_file.c_str()); for(int j=0; j<3; j++) { for(unsigned int i=0; iSize(); i++) { GradientDirectionType direction = gradients->ElementAt(i); direction.normalize(); myfile2 << direction.get(j) << " "; } myfile2 << std::endl; } } std::vector mitk::gradients::GetAllUniqueDirections(const BValueMap & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) { IndiciesVector directioncontainer; auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); auto containerIt = directioncontainer.begin(); bool directionExist = false; while(containerIt != directioncontainer.end()) { if (fabs(dot_product(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { directioncontainer.push_back(wntIndex); } } } return directioncontainer; } bool mitk::gradients::CheckForDifferingShellDirections(const BValueMap & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ auto mapIterator_2 = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator_2++; //skip bzero Values for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){ if(mapIterator_2 == mapIterator) continue; IndiciesVector currentShell = mapIterator->second; IndiciesVector testShell = mapIterator_2->second; for (unsigned int i = 0; i< currentShell.size(); i++) if (fabs(dot_product(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) { return true; } } } return false; } vnl_matrix mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer) { vnl_matrix Q(3, refShell.size()); Q.fill(0.0); for(unsigned int i = 0; i < refShell.size(); i++) { GradientDirectionType dir = refGradientsContainer->ElementAt(refShell[i]); double x = dir.normalize().get(0); double y = dir.normalize().get(1); double z = dir.normalize().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i) = cart[2]; } return Q; } vnl_matrix mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder) { vnl_matrix SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5); SHBasisOutput.fill(0.0); for(int i=0; i< (int)SHBasisOutput.rows(); i++) for(int k = 0; k <= (int)LOrder; k += 2) for(int m =- k; m <= k; m++) { int j = ( k * k + k + 2 ) / 2.0 + m - 1; double phi = QBallReference(0,i); double th = QBallReference(1,i); double val = mitk::sh::Yj(m,k,th,phi); SHBasisOutput(i,j) = val; } return SHBasisOutput; } mitk::gradients::GradientDirectionContainerType::Pointer mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(const mitk::gradients::BValueMap & bValueMap, const GradientDirectionContainerType *origninalGradentcontainer) { mitk::gradients::GradientDirectionContainerType::Pointer directioncontainer = mitk::gradients::GradientDirectionContainerType::New(); auto mapIterator = bValueMap.begin(); if(bValueMap.find(0) != bValueMap.end() && bValueMap.size() > 1){ mapIterator++; //skip bzero Values vnl_vector_fixed vec; vec.fill(0.0); directioncontainer->push_back(vec); } for( ; mapIterator != bValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin(); bool directionExist = false; while(containerIt != directioncontainer->End()) { if (fabs(dot_product(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex)); directioncontainer->push_back(dir.normalize()); } } } return directioncontainer; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp index dc28e08540..b2da19de97 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp @@ -1,446 +1,496 @@ /*=================================================================== 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 #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}) , m_SkipSelfConnections(false) , m_OnlySelfConnections(false) + , m_SplitByRoi(false) , m_SplitLabels(false) , m_MinFibersPerTract(0) , m_PairedStartEndLabels(false) { m_Interpolator = itk::LinearInterpolateImageFunction< itk::Image< PixelType, 3 >, float >::New(); } 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 > void FiberExtractionFilter< PixelType >::SetRoiImageNames(const std::vector< std::string > roi_names) { m_RoiImageNames = roi_names; } 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) { if( m_InputType == INPUT::SCALAR_MAP ) return mitk::imv::IsInsideMask(itkP, m_Interpolate, m_Interpolator, m_Threshold); else if( m_InputType == INPUT::LABEL_MAP ) { auto val = mitk::imv::GetImageValue(itkP, m_Interpolate, m_Interpolator); for (auto l : m_Labels) if (l==val) return true; } else mitkThrow() << "No valid input type selected!"; return false; } template< class PixelType > std::vector< std::string > FiberExtractionFilter< PixelType >::GetPositiveLabels() const { return m_PositiveLabels; } 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; + std::vector< std::vector< long > > positive_ids; // one ID vector per ROI 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; mSetInputImage(roi); - int inside = 0; - int outside = 0; - for (int j=0; jGetPoint(j); - - itk::Point itkP; - itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; - - if ( IsPositive(itkP) ) - inside++; - else - outside++; - - if ((float)inside/numPoints > m_OverlapFraction) + itk::Point startVertex = mitk::imv::GetItkPoint(points->GetPoint(j)); + itk::Index<3> startIndex; + itk::ContinuousIndex startIndexCont; + roi->TransformPhysicalPointToIndex(startVertex, startIndex); + roi->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); + + itk::Point endVertex = mitk::imv::GetItkPoint(points->GetPoint(j + 1)); + itk::Index<3> endIndex; + itk::ContinuousIndex endIndexCont; + roi->TransformPhysicalPointToIndex(endVertex, endIndex); + roi->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); + + std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(roi->GetSpacing(), startIndex, endIndex, startIndexCont, endIndexCont); + for (std::pair< itk::Index<3>, double > segment : segments) { - positive = true; - positive_ids[m].push_back(i); - break; + if (roi->GetPixel(segment.first)>=m_Threshold) + inside += segment.second; + else + outside += segment.second; } } + + if ((float)inside/(inside+outside) >= m_OverlapFraction) + { + 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) - if (ids.size()>=m_MinFibersPerTract) - m_Positives.push_back(CreateFib(ids)); + { + for (unsigned int i=0; i(i); + + if (positive_ids.at(i).size()>=m_MinFibersPerTract) + { + m_Positives.push_back(CreateFib(positive_ids.at(i))); + m_PositiveLabels.push_back(name); + } + } + if (!m_SplitByRoi) + { + mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); + output = output->AddBundles(m_Positives); + m_Positives.clear(); + m_Positives.push_back(output); + m_PositiveLabels.clear(); + m_PositiveLabels.push_back(""); + } + } } 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; + std::vector< std::vector< long > > positive_ids; // one ID vector per ROI 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; mSetInputImage(roi); int inside = 0; // check first fiber point { double* p = points->GetPoint(0); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; if ( IsPositive(itkP) ) 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]; if ( IsPositive(itkP) ) 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) - if (ids.size()>=m_MinFibersPerTract) - m_Positives.push_back(CreateFib(ids)); + { + for (unsigned int i=0; i(i); + + if (positive_ids.at(i).size()>=m_MinFibersPerTract) + { + m_Positives.push_back(CreateFib(positive_ids.at(i))); + m_PositiveLabels.push_back(name); + } + } + if (!m_SplitByRoi) + { + mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); + output = output->AddBundles(m_Positives); + m_Positives.clear(); + m_Positives.push_back(output); + m_PositiveLabels.clear(); + m_PositiveLabels.push_back(""); + } + } } template< class PixelType > void FiberExtractionFilter< PixelType >::ExtractLabels(mitk::FiberBundle::Pointer fib) { MITK_INFO << "Extracting fibers by labels"; vtkSmartPointer polydata = fib->GetFiberPolyData(); - std::vector< std::map< std::string, std::vector< long > > > positive_ids; - positive_ids.resize(m_RoiImages.size()); + std::map< std::string, std::vector< long > > positive_ids; std::vector< long > negative_ids; // fibers not overlapping with ANY label 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; mSetInputImage(roi); int inside = 0; double* p1 = points->GetPoint(0); itk::Point itkP1; itkP1[0] = p1[0]; itkP1[1] = p1[1]; itkP1[2] = p1[2]; short label1 = mitk::imv::GetImageValue(itkP1, m_Interpolate, m_Interpolator); double* p2 = points->GetPoint(numPoints-1); itk::Point itkP2; itkP2[0] = p2[0]; itkP2[1] = p2[1]; itkP2[2] = p2[2]; short label2 = mitk::imv::GetImageValue(itkP2, m_Interpolate, m_Interpolator); if (!m_Labels.empty()) // extract fibers from all pairwise label combinations { for (auto l : m_Labels) { if (l==label1) inside++; if (l==label2) inside++; if (inside==2) break; } } else // extract fibers between start and end labels { m_BothEnds = true; // if we have start and end labels it does not make sense to not use both endpoints if (m_PairedStartEndLabels) { if (m_StartLabels.size()!=m_EndLabels.size()) mitkThrow() << "Start and end label lists must have same size if paired labels are used"; for (unsigned int ii=0; ii(label1) + "-" + boost::lexical_cast(label2); - else - key = boost::lexical_cast(label2) + "-" + boost::lexical_cast(label1); - if (m(label1) + "-" + boost::lexical_cast(label2); + else + key = boost::lexical_cast(label2) + "-" + boost::lexical_cast(label1); + } + + if (m_SplitByRoi) + { + if (m(m) + "_" + key; + } if (m_BothEnds) { if ( (inside==2 && (!m_SkipSelfConnections || label1!=label2)) || (inside==2 && m_OnlySelfConnections && label1==label2) ) { positive = true; - if ( positive_ids[m].count(key)==0 ) - positive_ids[m].insert( std::pair< std::string, std::vector< long > >( key, {i}) ); + if ( positive_ids.count(key)==0 ) + positive_ids.insert( std::pair< std::string, std::vector< long > >( key, {i}) ); else - positive_ids[m][key].push_back(i); + positive_ids[key].push_back(i); } } else { if ( (inside>=1 && (!m_SkipSelfConnections || label1!=label2)) || (inside==2 && m_OnlySelfConnections && label1==label2) ) { positive = true; - if ( positive_ids[m].count(key)==0 ) - positive_ids[m].insert( std::pair< std::string, std::vector< long > >( key, {i}) ); + if ( positive_ids.count(key)==0 ) + positive_ids.insert( std::pair< std::string, std::vector< long > >( key, {i}) ); else - positive_ids[m][key].push_back(i); + positive_ids[key].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 labels : positive_ids) - for (auto tuple : labels) + for (auto label : positive_ids) + { + if (label.second.size()>=m_MinFibersPerTract) { - if (tuple.second.size()>=m_MinFibersPerTract) - { - m_Positives.push_back(CreateFib(tuple.second)); - m_PositiveLabels.push_back(tuple.first); - } + m_Positives.push_back(CreateFib(label.second)); + m_PositiveLabels.push_back(label.first); } - if (!m_SplitLabels) - { - mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); - output = output->AddBundles(m_Positives); - m_Positives.clear(); - m_Positives.push_back(output); - m_PositiveLabels.clear(); } + +// if (!m_SplitLabels) +// { +// mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); +// output = output->AddBundles(m_Positives); +// m_Positives.clear(); +// m_Positives.push_back(output); +// m_PositiveLabels.clear(); +// } } } template< class PixelType > void FiberExtractionFilter< PixelType >::SetLabels(const std::vector &Labels) { m_Labels = Labels; } template< class PixelType > void FiberExtractionFilter< PixelType >::SetStartLabels(const std::vector &Labels) { m_StartLabels = Labels; } template< class PixelType > void FiberExtractionFilter< PixelType >::SetEndLabels(const std::vector &Labels) { m_EndLabels = 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) { if (m_InputType==INPUT::LABEL_MAP) ExtractLabels(fib); else ExtractEndpoints(fib); } } } #endif // __itkFiberExtractionFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h index 24d2f1bbf2..da31262217 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.h @@ -1,131 +1,133 @@ /*=================================================================== 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 #include namespace itk{ /** * \brief Extract streamlines from tractograms using ROI images */ template< class PixelType > class FiberExtractionFilter : public ProcessObject { public: enum MODE { OVERLAP, ENDPOINTS }; enum INPUT { SCALAR_MAP, ///< In this case, positive means roi image vlaue > threshold LABEL_MAP ///< In this case, positive means roi image value 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 ) void Update() override{ this->GenerateData(); } itkSetMacro( InputFiberBundle, mitk::FiberBundle::Pointer ) itkSetMacro( Mode, MODE ) ///< Overlap or endpoints itkSetMacro( InputType, INPUT ) ///< Scalar map or label image itkSetMacro( BothEnds, bool ) ///< Both streamline ends need to be inside of the ROI for the streamline to be considered positive itkSetMacro( OverlapFraction, float ) ///< Necessary fraction of streamline points inside the ROI for the streamline to be considered positive itkSetMacro( DontResampleFibers, bool ) ///< Don't resample input fibers to ensure coverage itkSetMacro( NoNegatives, bool ) ///< Don't create output tractograms from negative streamlines (save the computation) itkSetMacro( NoPositives, bool ) ///< Don't create output tractograms from positive streamlines (save the computation) - itkSetMacro( Interpolate, bool ) ///< Interpolate input ROI image + itkSetMacro( Interpolate, bool ) ///< Interpolate input ROI image (only relevant for ENDPOINTS mode itkSetMacro( Threshold, float ) ///< Threshold on input ROI image value to determine positives/negatives itkSetMacro( SkipSelfConnections, bool ) ///< Ignore streamlines between two identical labels itkSetMacro( OnlySelfConnections, bool ) ///< Only keep streamlines between two identical labels itkSetMacro( SplitLabels, bool ) ///< Output a separate tractogram for each label-->label tract + itkSetMacro( SplitByRoi, bool ) ///< Output a separate tractogram for each ROI image itkSetMacro( MinFibersPerTract, unsigned int ) ///< Discard positives with less fibers itkSetMacro( PairedStartEndLabels, bool ) void SetRoiImages(const std::vector< ItkInputImgType* > &rois); void SetRoiImageNames(const std::vector< std::string > roi_names); void SetLabels(const std::vector &Labels); void SetStartLabels(const std::vector &Labels); void SetEndLabels(const std::vector &Labels); std::vector GetPositives() const; ///< Get positive tracts (filtered by the input ROIs) std::vector GetNegatives() const; ///< Get negative tracts (not filtered by the ROIs) std::vector< std::string > GetPositiveLabels() const; ///< In case of label extraction, this vector contains the labels corresponding to the positive tracts protected: void GenerateData() override; FiberExtractionFilter(); ~FiberExtractionFilter() override; mitk::FiberBundle::Pointer CreateFib(std::vector< long >& ids); void ExtractOverlap(mitk::FiberBundle::Pointer fib); void ExtractEndpoints(mitk::FiberBundle::Pointer fib); void ExtractLabels(mitk::FiberBundle::Pointer fib); bool IsPositive(const itk::Point& itkP); mitk::FiberBundle::Pointer m_InputFiberBundle; std::vector< mitk::FiberBundle::Pointer > m_Positives; std::vector< mitk::FiberBundle::Pointer > m_Negatives; std::vector< ItkInputImgType* > m_RoiImages; std::vector< std::string > m_RoiImageNames; 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; bool m_SkipSelfConnections; bool m_OnlySelfConnections; + bool m_SplitByRoi; bool m_SplitLabels; unsigned int m_MinFibersPerTract; std::vector< unsigned short > m_StartLabels; std::vector< unsigned short > m_EndLabels; bool m_PairedStartEndLabels; std::vector< std::string > m_PositiveLabels; typename itk::LinearInterpolateImageFunction< itk::Image< PixelType, 3 >, float >::Pointer m_Interpolator; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkFiberExtractionFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp index 3c7f599a2d..cc815686fe 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp @@ -1,997 +1,993 @@ #include "itkFitFibersToImageFilter.h" #include #include namespace itk{ FitFibersToImageFilter::FitFibersToImageFilter() : m_PeakImage(nullptr) , m_DiffImage(nullptr) , m_ScalarImage(nullptr) , m_MaskImage(nullptr) , m_FitIndividualFibers(true) , m_GradientTolerance(1e-5) , m_Lambda(0.1) , m_MaxIterations(20) , m_Coverage(0) , m_Overshoot(0) , m_RMSE(0.0) , m_FilterOutliers(false) , m_MeanWeight(1.0) , m_MedianWeight(1.0) , m_MinWeight(1.0) , m_MaxWeight(1.0) , m_Verbose(true) , m_NumUnknowns(0) , m_NumResiduals(0) , m_NumCoveredDirections(0) , m_SignalModel(nullptr) , sz_x(0) , sz_y(0) , sz_z(0) , m_MeanTractDensity(0) , m_MeanSignal(0) , fiber_count(0) , m_Regularization(VnlCostFunction::REGU::VOXEL_VARIANCE) { this->SetNumberOfRequiredOutputs(3); } FitFibersToImageFilter::~FitFibersToImageFilter() { } itk::Point FitFibersToImageFilter::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } void FitFibersToImageFilter::CreateDiffSystem() { sz_x = m_DiffImage->GetLargestPossibleRegion().GetSize(0); sz_y = m_DiffImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_DiffImage->GetLargestPossibleRegion().GetSize(2); dim_four_size = m_DiffImage->GetVectorLength(); int num_voxels = sz_x*sz_y*sz_z; auto spacing = m_DiffImage->GetSpacing(); m_NumResiduals = num_voxels * dim_four_size; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; m_NumCoveredDirections = 0; fiber_count = 0; vnl_vector voxel_indicator; voxel_indicator.set_size(sz_x*sz_y*sz_z); voxel_indicator.fill(0); int numFibers = 0; for (unsigned int bundle=0; bundleGetNumFibers(); boost::progress_display disp(numFibers); m_GroupSizes.clear(); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers()); for (int i=0; iGetNumFibers(); ++i) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_DiffImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_DiffImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); PointType3 endVertex = GetItkPoint(points->GetPoint(j+1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_DiffImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_DiffImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); mitk::DiffusionSignalModel<>::GradientType fiber_dir; fiber_dir[0] = endVertex[0]-startVertex[0]; fiber_dir[1] = endVertex[1]-startVertex[1]; fiber_dir[2] = endVertex[2]-startVertex[2]; fiber_dir.Normalize(); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_DiffImage->GetLargestPossibleRegion().IsInside(seg.first) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(seg.first)==0)) continue; int x = seg.first[0]; int y = seg.first[1]; int z = seg.first[2]; mitk::DiffusionSignalModel<>::PixelType simulated_pixel = m_SignalModel->SimulateMeasurement(fiber_dir)*seg.second; VectorImgType::PixelType measured_pixel = m_DiffImage->GetPixel(seg.first); double simulated_mean = 0; double measured_mean = 0; int num_nonzero_g = 0; for (int g=0; gGetGradientDirection(g).GetNorm()GetLargestPossibleRegion().GetSize(0); sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2); dim_four_size = m_PeakImage->GetLargestPossibleRegion().GetSize(3)/3 + 1; // +1 for zero - peak int num_voxels = sz_x*sz_y*sz_z; itk::Vector< double, 3 > spacing; spacing[0] = m_PeakImage->GetSpacing()[0]; spacing[1] = m_PeakImage->GetSpacing()[1]; spacing[2] = m_PeakImage->GetSpacing()[2]; PointType3 origin; origin[0] = m_PeakImage->GetOrigin()[0]; origin[1] = m_PeakImage->GetOrigin()[1]; origin[2] = m_PeakImage->GetOrigin()[2]; UcharImgType::RegionType::SizeType size; size[0] = m_PeakImage->GetLargestPossibleRegion().GetSize()[0]; size[1] = m_PeakImage->GetLargestPossibleRegion().GetSize()[1]; size[2] = m_PeakImage->GetLargestPossibleRegion().GetSize()[2]; UcharImgType::RegionType imageRegion; imageRegion.SetSize(size); itk::Matrix direction; for (int r=0; r<3; r++) for (int c=0; c<3; c++) direction[r][c] = m_PeakImage->GetDirection()[r][c]; if (m_MaskImage.IsNull()) { m_MaskImage = UcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_NumResiduals = num_voxels * dim_four_size; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; m_NumCoveredDirections = 0; fiber_count = 0; m_GroupSizes.clear(); int numFibers = 0; for (unsigned int bundle=0; bundleGetNumFibers(); boost::progress_display disp(numFibers); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers()); for (int i=0; iGetNumFibers(); ++i) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_MaskImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_MaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); PointType3 endVertex = GetItkPoint(points->GetPoint(j+1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_MaskImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_MaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); vnl_vector_fixed fiber_dir; fiber_dir[0] = endVertex[0]-startVertex[0]; fiber_dir[1] = endVertex[1]-startVertex[1]; fiber_dir[2] = endVertex[2]-startVertex[2]; fiber_dir.normalize(); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_MaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_MaskImage->GetPixel(seg.first)==0) continue; itk::Index<4> idx4; idx4[0]=seg.first[0]; idx4[1]=seg.first[1]; idx4[2]=seg.first[2]; idx4[3]=0; double w = 1; int peak_id = dim_four_size-1; double peak_mag = 0; GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w, peak_mag); w *= seg.second; int x = idx4[0]; int y = idx4[1]; int z = idx4[2]; unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*peak_id; if (b[linear_index] == 0 && peak_idGetLargestPossibleRegion().GetSize(0); sz_y = m_ScalarImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_ScalarImage->GetLargestPossibleRegion().GetSize(2); int num_voxels = sz_x*sz_y*sz_z; auto spacing = m_ScalarImage->GetSpacing(); m_NumResiduals = num_voxels; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; int numCoveredVoxels = 0; fiber_count = 0; int numFibers = 0; for (unsigned int bundle=0; bundleGetNumFibers(); boost::progress_display disp(numFibers); m_GroupSizes.clear(); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); m_GroupSizes.push_back(m_Tractograms.at(bundle)->GetNumFibers()); for (int i=0; iGetNumFibers(); ++i) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_ScalarImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_ScalarImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); PointType3 endVertex = GetItkPoint(points->GetPoint(j+1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_ScalarImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_ScalarImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_ScalarImage->GetLargestPossibleRegion().IsInside(seg.first) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(seg.first)==0)) continue; float image_value = m_ScalarImage->GetPixel(seg.first); int x = seg.first[0]; int y = seg.first[1]; int z = seg.first[2]; unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z; if (b[linear_index] == 0) { numCoveredVoxels++; m_MeanSignal += image_value; } m_MeanTractDensity += seg.second; if (m_FitIndividualFibers) { b[linear_index] = image_value; A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + seg.second); } else { b[linear_index] = image_value; A.put(linear_index, bundle, A.get(linear_index, bundle) + seg.second); } } } ++fiber_count; } } m_MeanTractDensity /= (numCoveredVoxels*fiber_count); m_MeanSignal /= numCoveredVoxels; A /= m_MeanTractDensity; b *= 100.0/m_MeanSignal; // times 100 because we want to avoid too small values for computational reasons // NEW FIT // m_MeanTractDensity /= numCoveredVoxels; // m_MeanSignal /= numCoveredVoxels; // b /= m_MeanSignal; // b *= m_MeanTractDensity; } void FitFibersToImageFilter::GenerateData() { m_NumUnknowns = m_Tractograms.size(); if (m_FitIndividualFibers) { m_NumUnknowns = 0; for (unsigned int bundle=0; bundleGetNumFibers(); } else m_FilterOutliers = false; if (m_NumUnknowns<1) { MITK_INFO << "No fibers in tractogram."; return; } fiber_count = 0; sz_x = 0; sz_y = 0; sz_z = 0; m_MeanTractDensity = 0; m_MeanSignal = 0; if (m_PeakImage.IsNotNull()) CreatePeakSystem(); else if (m_DiffImage.IsNotNull()) CreateDiffSystem(); else if (m_ScalarImage.IsNotNull()) CreateScalarSystem(); else mitkThrow() << "No input image set!"; MITK_INFO << "Initializing optimizer"; double init_lambda = fiber_count; // initialization for lambda estimation itk::TimeProbe clock; clock.Start(); cost = VnlCostFunction(m_NumUnknowns); cost.SetProblem(A, b, init_lambda, m_Regularization); cost.SetGroupSizes(m_GroupSizes); m_Weights.set_size(m_NumUnknowns); m_Weights.fill( 1.0/m_NumUnknowns ); vnl_lbfgsb minimizer(cost); vnl_vector l; l.set_size(m_NumUnknowns); l.fill(0); vnl_vector bound_selection; bound_selection.set_size(m_NumUnknowns); bound_selection.fill(1); minimizer.set_bound_selection(bound_selection); minimizer.set_lower_bound(l); minimizer.set_projected_gradient_tolerance(m_GradientTolerance); if (m_Regularization==VnlCostFunction::REGU::MSM) MITK_INFO << "Regularization type: MSM"; else if (m_Regularization==VnlCostFunction::REGU::VARIANCE) MITK_INFO << "Regularization type: VARIANCE"; else if (m_Regularization==VnlCostFunction::REGU::LASSO) MITK_INFO << "Regularization type: LASSO"; else if (m_Regularization==VnlCostFunction::REGU::VOXEL_VARIANCE) MITK_INFO << "Regularization type: VOXEL_VARIANCE"; else if (m_Regularization==VnlCostFunction::REGU::GROUP_LASSO) MITK_INFO << "Regularization type: GROUP_LASSO"; else if (m_Regularization==VnlCostFunction::REGU::GROUP_VARIANCE) MITK_INFO << "Regularization type: GROUP_VARIANCE"; else if (m_Regularization==VnlCostFunction::REGU::NONE) MITK_INFO << "Regularization type: NONE"; if (m_Regularization!=VnlCostFunction::REGU::NONE) // REMOVE FOR NEW FIT AND SET cost.m_Lambda = m_Lambda { MITK_INFO << "Estimating regularization"; minimizer.set_trace(false); minimizer.set_max_function_evals(2); minimizer.minimize(m_Weights); vnl_vector dx; dx.set_size(m_NumUnknowns); dx.fill(0.0); cost.calc_regularization_gradient(m_Weights, dx); if (m_Weights.magnitude()==0) { MITK_INFO << "Regularization estimation failed. Using default value."; cost.m_Lambda = fiber_count*m_Lambda; } else { double r = dx.magnitude()/m_Weights.magnitude(); // wtf??? cost.m_Lambda *= m_Lambda*55.0/r; MITK_INFO << r << " - " << m_Lambda*55.0/r; if (cost.m_Lambda>10e7) { MITK_INFO << "Regularization estimation failed. Using default value."; cost.m_Lambda = fiber_count*m_Lambda; } } } else cost.m_Lambda = 0; MITK_INFO << "Using regularization factor of " << cost.m_Lambda << " (λ: " << m_Lambda << ")"; MITK_INFO << "Fitting fibers"; minimizer.set_trace(m_Verbose); minimizer.set_max_function_evals(m_MaxIterations); minimizer.minimize(m_Weights); std::vector< double > weights; if (m_FilterOutliers) { for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); MITK_INFO << "Setting upper weight bound to " << weights.at(m_NumUnknowns*0.99); vnl_vector u; u.set_size(m_NumUnknowns); u.fill(weights.at(m_NumUnknowns*0.99)); minimizer.set_upper_bound(u); bound_selection.fill(2); minimizer.set_bound_selection(bound_selection); minimizer.minimize(m_Weights); weights.clear(); } for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); m_MeanWeight = m_Weights.mean(); m_MedianWeight = weights.at(m_NumUnknowns*0.5); m_MinWeight = weights.at(0); m_MaxWeight = weights.at(m_NumUnknowns-1); MITK_INFO << "*************************"; MITK_INFO << "Weight statistics"; MITK_INFO << "Sum: " << m_Weights.sum(); MITK_INFO << "Mean: " << m_MeanWeight; MITK_INFO << "1% quantile: " << weights.at(m_NumUnknowns*0.01); MITK_INFO << "5% quantile: " << weights.at(m_NumUnknowns*0.05); MITK_INFO << "25% quantile: " << weights.at(m_NumUnknowns*0.25); MITK_INFO << "Median: " << m_MedianWeight; MITK_INFO << "75% quantile: " << weights.at(m_NumUnknowns*0.75); MITK_INFO << "95% quantile: " << weights.at(m_NumUnknowns*0.95); MITK_INFO << "99% quantile: " << weights.at(m_NumUnknowns*0.99); MITK_INFO << "Min: " << m_MinWeight; MITK_INFO << "Max: " << m_MaxWeight; MITK_INFO << "*************************"; MITK_INFO << "NumEvals: " << minimizer.get_num_evaluations(); MITK_INFO << "NumIterations: " << minimizer.get_num_iterations(); MITK_INFO << "Residual cost: " << minimizer.get_end_error(); m_RMSE = cost.S->get_rms_error(m_Weights); MITK_INFO << "Final RMSE: " << m_RMSE; clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Optimization took " << h << "h, " << m << "m and " << s << "s"; MITK_INFO << "Weighting fibers"; m_RmsDiffPerBundle.set_size(m_Tractograms.size()); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); if (m_FitIndividualFibers) { unsigned int fiber_count = 0; for (unsigned int bundle=0; bundle temp_weights; temp_weights.set_size(m_Weights.size()); temp_weights.copy_in(m_Weights.data_block()); for (int i=0; iGetNumFibers(); i++) { m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]); temp_weights[fiber_count] = 0; ++fiber_count; } double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE; m_RmsDiffPerBundle[bundle] = d_rms; - m_Tractograms.at(bundle)->Compress(0.1); - m_Tractograms.at(bundle)->ColorFibersByFiberWeights(false, true); } } else { for (unsigned int i=0; i temp_weights; temp_weights.set_size(m_Weights.size()); temp_weights.copy_in(m_Weights.data_block()); temp_weights[i] = 0; double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE; m_RmsDiffPerBundle[i] = d_rms; m_Tractograms.at(i)->SetFiberWeights(m_Weights[i]); - m_Tractograms.at(i)->Compress(0.1); - m_Tractograms.at(i)->ColorFibersByFiberWeights(false, true); } } std::cout.rdbuf (old); // transform back A *= m_MeanSignal/100.0; b *= m_MeanSignal/100.0; MITK_INFO << "Generating output images ..."; if (m_PeakImage.IsNotNull()) GenerateOutputPeakImages(); else if (m_DiffImage.IsNotNull()) GenerateOutputDiffImages(); else if (m_ScalarImage.IsNotNull()) GenerateOutputScalarImages(); m_Coverage = m_Coverage/m_MeanSignal; m_Overshoot = m_Overshoot/m_MeanSignal; MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*m_Coverage << "%"; MITK_INFO << std::fixed << "Overshoot: " << setprecision(2) << 100.0*m_Overshoot << "%"; } void FitFibersToImageFilter::GenerateOutputDiffImages() { VectorImgType::PixelType pix; pix.SetSize(m_DiffImage->GetVectorLength()); pix.Fill(0); itk::ImageDuplicator< VectorImgType >::Pointer duplicator = itk::ImageDuplicator< VectorImgType >::New(); duplicator->SetInputImage(m_DiffImage); duplicator->Update(); m_UnderexplainedImageDiff = duplicator->GetOutput(); m_UnderexplainedImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_UnderexplainedImageDiff); duplicator->Update(); m_OverexplainedImageDiff = duplicator->GetOutput(); m_OverexplainedImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_OverexplainedImageDiff); duplicator->Update(); m_ResidualImageDiff = duplicator->GetOutput(); m_ResidualImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_ResidualImageDiff); duplicator->Update(); m_FittedImageDiff = duplicator->GetOutput(); m_FittedImageDiff->FillBuffer(pix); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); itk::ImageRegionIterator it1 = itk::ImageRegionIterator(m_DiffImage, m_DiffImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2 = itk::ImageRegionIterator(m_FittedImageDiff, m_FittedImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it3 = itk::ImageRegionIterator(m_ResidualImageDiff, m_ResidualImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it4 = itk::ImageRegionIterator(m_UnderexplainedImageDiff, m_UnderexplainedImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it5 = itk::ImageRegionIterator(m_OverexplainedImageDiff, m_OverexplainedImageDiff->GetLargestPossibleRegion()); m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; while( !it2.IsAtEnd() ) { itk::Index<3> idx3 = it2.GetIndex(); VectorImgType::PixelType original_pix =it1.Get(); VectorImgType::PixelType fitted_pix =it2.Get(); VectorImgType::PixelType residual_pix =it3.Get(); VectorImgType::PixelType underexplained_pix =it4.Get(); VectorImgType::PixelType overexplained_pix =it5.Get(); int num_nonzero_g = 0; double original_mean = 0; for (int g=0; gGetGradientDirection(g).GetNorm()>=mitk::eps ) { original_mean += original_pix[g]; ++num_nonzero_g; } } original_mean /= num_nonzero_g; for (int g=0; g=0) { underexplained_pix[g] = residual_pix[g]; m_Coverage += fitted_b[linear_index] + original_mean; } m_MeanSignal += b[linear_index] + original_mean; } it2.Set(fitted_pix); it3.Set(residual_pix); it4.Set(underexplained_pix); it5.Set(overexplained_pix); ++it1; ++it2; ++it3; ++it4; ++it5; } } void FitFibersToImageFilter::GenerateOutputScalarImages() { itk::ImageDuplicator< DoubleImgType >::Pointer duplicator = itk::ImageDuplicator< DoubleImgType >::New(); duplicator->SetInputImage(m_ScalarImage); duplicator->Update(); m_UnderexplainedImageScalar = duplicator->GetOutput(); m_UnderexplainedImageScalar->FillBuffer(0); duplicator->SetInputImage(m_UnderexplainedImageScalar); duplicator->Update(); m_OverexplainedImageScalar = duplicator->GetOutput(); m_OverexplainedImageScalar->FillBuffer(0); duplicator->SetInputImage(m_OverexplainedImageScalar); duplicator->Update(); m_ResidualImageScalar = duplicator->GetOutput(); m_ResidualImageScalar->FillBuffer(0); duplicator->SetInputImage(m_ResidualImageScalar); duplicator->Update(); m_FittedImageScalar = duplicator->GetOutput(); m_FittedImageScalar->FillBuffer(0); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); itk::ImageRegionIterator it1 = itk::ImageRegionIterator(m_ScalarImage, m_ScalarImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2 = itk::ImageRegionIterator(m_FittedImageScalar, m_FittedImageScalar->GetLargestPossibleRegion()); itk::ImageRegionIterator it3 = itk::ImageRegionIterator(m_ResidualImageScalar, m_ResidualImageScalar->GetLargestPossibleRegion()); itk::ImageRegionIterator it4 = itk::ImageRegionIterator(m_UnderexplainedImageScalar, m_UnderexplainedImageScalar->GetLargestPossibleRegion()); itk::ImageRegionIterator it5 = itk::ImageRegionIterator(m_OverexplainedImageScalar, m_OverexplainedImageScalar->GetLargestPossibleRegion()); m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; while( !it2.IsAtEnd() ) { itk::Index<3> idx3 = it2.GetIndex(); DoubleImgType::PixelType original_pix =it1.Get(); DoubleImgType::PixelType fitted_pix =it2.Get(); DoubleImgType::PixelType residual_pix =it3.Get(); DoubleImgType::PixelType underexplained_pix =it4.Get(); DoubleImgType::PixelType overexplained_pix =it5.Get(); unsigned int linear_index = idx3[0] + sz_x*idx3[1] + sz_x*sz_y*idx3[2]; fitted_pix = fitted_b[linear_index]; residual_pix = original_pix - fitted_pix; if (residual_pix<0) { overexplained_pix = residual_pix; m_Coverage += b[linear_index]; m_Overshoot -= residual_pix; } else if (residual_pix>=0) { underexplained_pix = residual_pix; m_Coverage += fitted_b[linear_index]; } m_MeanSignal += b[linear_index]; it2.Set(fitted_pix); it3.Set(residual_pix); it4.Set(underexplained_pix); it5.Set(overexplained_pix); ++it1; ++it2; ++it3; ++it4; ++it5; } } VnlCostFunction::REGU FitFibersToImageFilter::GetRegularization() const { return m_Regularization; } void FitFibersToImageFilter::SetRegularization(const VnlCostFunction::REGU &Regularization) { m_Regularization = Regularization; } void FitFibersToImageFilter::GenerateOutputPeakImages() { itk::ImageDuplicator< PeakImgType >::Pointer duplicator = itk::ImageDuplicator< PeakImgType >::New(); duplicator->SetInputImage(m_PeakImage); duplicator->Update(); m_UnderexplainedImage = duplicator->GetOutput(); m_UnderexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_UnderexplainedImage); duplicator->Update(); m_OverexplainedImage = duplicator->GetOutput(); m_OverexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_OverexplainedImage); duplicator->Update(); m_ResidualImage = duplicator->GetOutput(); m_ResidualImage->FillBuffer(0.0); duplicator->SetInputImage(m_ResidualImage); duplicator->Update(); m_FittedImage = duplicator->GetOutput(); m_FittedImage->FillBuffer(0.0); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); for (unsigned int r=0; r idx4; unsigned int linear_index = r; idx4[0] = linear_index % sz_x; linear_index /= sz_x; idx4[1] = linear_index % sz_y; linear_index /= sz_y; idx4[2] = linear_index % sz_z; linear_index /= sz_z; int peak_id = linear_index % dim_four_size; if (peak_id peak_dir; idx4[3] = peak_id*3; peak_dir[0] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[1] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[2] = m_PeakImage->GetPixel(idx4); peak_dir.normalize(); peak_dir *= fitted_b[r]; idx4[3] = peak_id*3; m_FittedImage->SetPixel(idx4, peak_dir[0]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[1]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[2]); } } m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; itk::Index<4> idx4; for (idx4[0]=0; idx4[0] idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2]; if (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0) continue; vnl_vector_fixed peak_dir; vnl_vector_fixed fitted_dir; vnl_vector_fixed overshoot_dir; for (idx4[3]=0; idx4[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx4[3]) { peak_dir[idx4[3]%3] = m_PeakImage->GetPixel(idx4); fitted_dir[idx4[3]%3] = m_FittedImage->GetPixel(idx4); m_ResidualImage->SetPixel(idx4, m_PeakImage->GetPixel(idx4) - m_FittedImage->GetPixel(idx4)); if (idx4[3]%3==2) { m_MeanSignal += peak_dir.magnitude(); itk::Index<4> tidx= idx4; if (peak_dir.magnitude()>fitted_dir.magnitude()) { m_Coverage += fitted_dir.magnitude(); m_UnderexplainedImage->SetPixel(tidx, peak_dir[2]-fitted_dir[2]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[1]-fitted_dir[1]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[0]-fitted_dir[0]); } else { overshoot_dir[0] = fitted_dir[0]-peak_dir[0]; overshoot_dir[1] = fitted_dir[1]-peak_dir[1]; overshoot_dir[2] = fitted_dir[2]-peak_dir[2]; m_Coverage += peak_dir.magnitude(); m_Overshoot += overshoot_dir.magnitude(); m_OverexplainedImage->SetPixel(tidx, overshoot_dir[2]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[1]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[0]); } } } } } void FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed fiber_dir, int& id, double& w, double& peak_mag ) { int m_NumDirs = peak_image->GetLargestPossibleRegion().GetSize()[3]/3; vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0.9; for (int i=0; i dir; idx[3] = i*3; dir[0] = peak_image->GetPixel(idx); idx[3] += 1; dir[1] = peak_image->GetPixel(idx); idx[3] += 1; dir[2] = peak_image->GetPixel(idx); float mag = dir.magnitude(); if (magangle) { angle = a; w = angle; peak_mag = mag; id = i; } } } std::vector FitFibersToImageFilter::GetTractograms() const { return m_Tractograms; } void FitFibersToImageFilter::SetTractograms(const std::vector &tractograms) { m_Tractograms = tractograms; } void FitFibersToImageFilter::SetSignalModel(mitk::DiffusionSignalModel<> *SignalModel) { m_SignalModel = SignalModel; } } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp index 86c7ec8a3f..00c8f7ac71 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp @@ -1,265 +1,261 @@ /*=================================================================== 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 ItkFloatImgType; ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = mitk::IOUtil::Load(filename); 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.beginGroup("1. Mandatory arguments:"); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false); parser.addArgument("rois", "", mitkCommandLineParser::StringList, "ROI images:", "ROI images", us::Any(), false); parser.endGroup(); parser.beginGroup("2. Label based extraction:"); parser.addArgument("labels", "", mitkCommandLineParser::StringList, "Labels:", "positive means roi image value in labels vector", us::Any()); parser.addArgument("split_labels", "", mitkCommandLineParser::Bool, "Split labels:", "output a separate tractogram for each label-->label tract", false); parser.addArgument("skip_self_connections", "", mitkCommandLineParser::Bool, "Skip self connections:", "ignore streamlines between two identical labels", false); parser.addArgument("start_labels", "", mitkCommandLineParser::StringList, "Start Labels:", "use separate start and end labels instead of one mixed set", us::Any()); parser.addArgument("end_labels", "", mitkCommandLineParser::StringList, "End Labels:", "use separate start and end labels instead of one mixed set", us::Any()); parser.addArgument("paired", "", mitkCommandLineParser::Bool, "Paired:", "start and end label list are paired", false); parser.endGroup(); parser.beginGroup("3. Misc:"); 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 ROI images", false); parser.addArgument("output_negatives", "", mitkCommandLineParser::Bool, "Negatives:", "output negatives", false); - parser.addArgument("interpolate", "", mitkCommandLineParser::Bool, "Interpolate:", "interpolate ROI images", false); + parser.addArgument("interpolate", "", mitkCommandLineParser::Bool, "Interpolate:", "interpolate ROI images (only for endpoint based extraction)", false); parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", 0.5); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. num. fibers:", "discard positive tracts with less fibers", 0); + parser.addArgument("split_rois", "", mitkCommandLineParser::Bool, "Split ROIs:", "output a separate tractogram for each ROI", false); parser.endGroup(); 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"]); unsigned int min_fibers = 0; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); bool split_labels = false; if (parsedArgs.count("split_labels")) split_labels = us::any_cast(parsedArgs["split_labels"]); + bool split_rois = false; + if (parsedArgs.count("split_rois")) + split_rois = us::any_cast(parsedArgs["split_rois"]); + bool skip_self_connections = false; if (parsedArgs.count("skip_self_connections")) skip_self_connections = us::any_cast(parsedArgs["skip_self_connections"]); bool output_negatives = false; if (parsedArgs.count("output_negatives")) output_negatives = us::any_cast(parsedArgs["output_negatives"]); 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"]); mitkCommandLineParser::StringContainerType labels; if (parsedArgs.count("labels")) labels = us::any_cast(parsedArgs["labels"]); mitkCommandLineParser::StringContainerType start_labels; if (parsedArgs.count("start_labels")) start_labels = us::any_cast(parsedArgs["start_labels"]); mitkCommandLineParser::StringContainerType end_labels; if (parsedArgs.count("end_labels")) end_labels = us::any_cast(parsedArgs["end_labels"]); bool paired = false; if (parsedArgs.count("paired")) paired = us::any_cast(parsedArgs["paired"]); try { // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(inFib); std::vector< ItkFloatImgType::Pointer > roi_images; std::vector< std::string > roi_names; for (std::size_t i=0; i roi_images2; for (auto roi : roi_images) roi_images2.push_back(roi); std::vector< unsigned short > short_labels; for (auto l : labels) short_labels.push_back(boost::lexical_cast(l)); std::vector< unsigned short > short_start_labels; for (auto l : start_labels) short_start_labels.push_back(boost::lexical_cast(l)); std::vector< unsigned short > short_end_labels; for (auto l : end_labels) short_end_labels.push_back(boost::lexical_cast(l)); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(inputTractogram); extractor->SetRoiImages(roi_images2); extractor->SetRoiImageNames(roi_names); extractor->SetOverlapFraction(overlap_fraction); extractor->SetBothEnds(both_ends); extractor->SetInterpolate(interpolate); extractor->SetThreshold(threshold); extractor->SetLabels(short_labels); extractor->SetStartLabels(short_start_labels); extractor->SetEndLabels(short_end_labels); extractor->SetSplitLabels(split_labels); + extractor->SetSplitByRoi(split_rois); extractor->SetMinFibersPerTract(min_fibers); extractor->SetSkipSelfConnections(skip_self_connections); extractor->SetPairedStartEndLabels(paired); if (!any_point) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); if (short_labels.size()>0 || short_start_labels.size()>0 || short_end_labels.size()>0) extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); extractor->Update(); - mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(nullptr); if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib + ".trk"); else { - if (!split_labels) - { - newFib = newFib->AddBundles(extractor->GetPositives()); - mitk::IOUtil::Save(newFib, outFib + ".trk"); - } - else + int c = 0; + std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); + for (auto fib : extractor->GetPositives()) { - int c = 0; - std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); - for (auto fib : extractor->GetPositives()) - { - std::string l = positive_labels.at(c); + std::string l = positive_labels.at(c); + if (l.size()>0) mitk::IOUtil::Save(fib, outFib + "_" + l + ".trk"); - ++c; - } + else + mitk::IOUtil::Save(fib, outFib + ".trk"); + ++c; } } if (output_negatives) { invert = !invert; if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib + "_negatives.trk"); else { - if (!split_labels) - { - newFib = newFib->AddBundles(extractor->GetPositives()); - mitk::IOUtil::Save(newFib, outFib + "_negatives.trk"); - } - else + int c = 0; + std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); + for (auto fib : extractor->GetPositives()) { - int c = 0; - std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); - for (auto fib : extractor->GetPositives()) - { - std::string l = positive_labels.at(c); + std::string l = positive_labels.at(c); + + if (l.size()>0) mitk::IOUtil::Save(fib, outFib + "_" + l + "_negatives.trk"); - ++c; - } + else + mitk::IOUtil::Save(fib, outFib + "_negatives.trk"); + ++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/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp index 0a616ad8bf..04f5e3e513 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp @@ -1,558 +1,594 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; std::vector< mitk::FiberBundle::Pointer > CombineTractograms(std::vector< mitk::FiberBundle::Pointer > reference, std::vector< mitk::FiberBundle::Pointer > candidates, int skip=-1) { std::vector< mitk::FiberBundle::Pointer > fib; for (auto f : reference) fib.push_back(f); int c = 0; for (auto f : candidates) { if (c!=skip) fib.push_back(f); ++c; } return fib; } 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 Score input candidate tracts using ACP analysis */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Anchor Constrained Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription("Score input candidate tracts using ACP analysis"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "a", mitkCommandLineParser::InputFile, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any(), false); parser.addArgument("", "p", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false); parser.addArgument("", "c", mitkCommandLineParser::InputDirectory, "Candidates folder:", "folder containing candidate tracts", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output folder:", "output folder", us::Any(), false); - parser.addArgument("anchor_masks", "", mitkCommandLineParser::StringList, "Reference Masks:", "reference tract masks for accuracy evaluation"); + parser.addArgument("reference_mask_folders", "", mitkCommandLineParser::StringList, "Reference Mask Folder(s):", "Folder(s) containing reference tract masks for accuracy evaluation"); parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "Mask image:", "scoring is only performed inside the mask image"); parser.addArgument("greedy_add", "", mitkCommandLineParser::Bool, "Greedy:", "if enabled, the candidate tracts are not jointly fitted to the residual image but one after the other employing a greedy scheme", false); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); parser.addArgument("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("regu", "", mitkCommandLineParser::String, "Regularization:", "MSM, Variance, VoxelVariance, Lasso, GroupLasso, GroupVariance, NONE (default)"); parser.addArgument("use_num_streamlines", "", mitkCommandLineParser::Bool, "Use number of streamlines as score:", "Don't fit candidates, simply use number of streamlines per candidate as score", false); parser.addArgument("use_weights", "", mitkCommandLineParser::Bool, "Use input weights as score:", "Don't fit candidates, simply use first input streamline weight per candidate as score", false); parser.addArgument("filter_zero_weights", "", mitkCommandLineParser::Bool, "Filter zero-weights", "Remove streamlines with weight 0 from candidates", false); + parser.addArgument("flipx", "", mitkCommandLineParser::Bool, "Flip x", "flip along x-axis", false); + parser.addArgument("flipy", "", mitkCommandLineParser::Bool, "Flip y", "flip along y-axis", false); + parser.addArgument("flipz", "", mitkCommandLineParser::Bool, "Flip z", "flip along z-axis", false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string anchors_file = us::any_cast(parsedArgs["a"]); std::string peak_file_name = us::any_cast(parsedArgs["p"]); std::string candidate_tract_folder = us::any_cast(parsedArgs["c"]); std::string out_folder = us::any_cast(parsedArgs["o"]); + if (!out_folder.empty() && out_folder.back() != '/') + out_folder += "/"; + if (!candidate_tract_folder.empty() && candidate_tract_folder.back() != '/') + candidate_tract_folder += "/"; + bool greedy_add = false; if (parsedArgs.count("greedy_add")) greedy_add = us::any_cast(parsedArgs["greedy_add"]); float lambda = 0.1; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); bool filter_zero_weights = false; if (parsedArgs.count("filter_zero_weights")) filter_zero_weights = us::any_cast(parsedArgs["filter_zero_weights"]); std::string mask_file = ""; if (parsedArgs.count("mask")) mask_file = us::any_cast(parsedArgs["mask"]); - mitkCommandLineParser::StringContainerType anchor_mask_files_folders; - if (parsedArgs.count("anchor_masks")) - anchor_mask_files_folders = us::any_cast(parsedArgs["anchor_masks"]); + mitkCommandLineParser::StringContainerType reference_mask_files_folders; + if (parsedArgs.count("reference_mask_folders")) + reference_mask_files_folders = us::any_cast(parsedArgs["reference_mask_folders"]); std::string regu = "NONE"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool use_weights = false; if (parsedArgs.count("use_weights")) use_weights = us::any_cast(parsedArgs["use_weights"]); bool use_num_streamlines = false; if (parsedArgs.count("use_num_streamlines")) use_num_streamlines = us::any_cast(parsedArgs["use_num_streamlines"]); + + bool flipx = false; + if (parsedArgs.count("flipx")) + flipx = us::any_cast(parsedArgs["flipx"]); + + bool flipy = false; + if (parsedArgs.count("flipy")) + flipy = us::any_cast(parsedArgs["flipy"]); + + bool flipz = false; + if (parsedArgs.count("z")) + flipz = us::any_cast(parsedArgs["flipz"]); + try { itk::TimeProbe clock; clock.Start(); if (!ist::PathExists(out_folder)) { MITK_INFO << "Creating output directory"; ist::MakeDirectory(out_folder); } MITK_INFO << "Loading data"; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect ofstream logfile; logfile.open (out_folder + "log.txt"); itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {}); mitk::Image::Pointer inputImage = dynamic_cast(mitk::IOUtil::Load(peak_file_name, &functor)[0].GetPointer()); // Load mask file. Fit is only performed inside the mask itk::FitFibersToImageFilter::UcharImgType::Pointer mask = nullptr; if (mask_file.compare("")!=0) { mitk::Image::Pointer mitk_mask = mitk::IOUtil::Load(mask_file); mitk::CastToItkImage(mitk_mask, mask); } // Load masks covering the true positives for evaluation purposes std::vector< itk::FitFibersToImageFilter::UcharImgType::Pointer > reference_masks; std::vector< std::string > anchor_mask_files; - for (auto filename : anchor_mask_files_folders) + for (auto filename : reference_mask_files_folders) { if (itksys::SystemTools::PathExists(filename)) { + if (!filename.empty() && filename.back() != '/') + filename += "/"; + auto list = get_file_list(filename, {".nrrd",".nii.gz",".nii"}); for (auto f : list) { MITK_INFO << f; itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr; mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load(f); mitk::CastToItkImage(ref_mitk_mask, ref_mask); reference_masks.push_back(ref_mask); anchor_mask_files.push_back(f); } } else if (itksys::SystemTools::FileExists(filename)) { anchor_mask_files.push_back(filename); itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr; mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load(filename); mitk::CastToItkImage(ref_mitk_mask, ref_mask); reference_masks.push_back(ref_mask); } } // Load peak image typedef mitk::ImageToItk< PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(inputImage); caster->Update(); PeakImgType::Pointer peak_image = caster->GetOutput(); // Load all candidate tracts std::vector< std::string > candidate_tract_files = get_file_list(candidate_tract_folder); std::vector< mitk::FiberBundle::Pointer > input_candidates; for (std::string f : candidate_tract_files) { mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(f); if (fib.IsNull()) continue; if (fib->GetNumFibers()<=0) continue; input_candidates.push_back(fib); } std::cout.rdbuf (old); // <-- restore MITK_INFO << "Loaded " << candidate_tract_files.size() << " candidate tracts."; MITK_INFO << "Loaded " << reference_masks.size() << " reference masks."; + if (flipx || flipy || flipz) + { + itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); + flipper->SetInput(peak_image); + flipper->SetFlipX(flipx); + flipper->SetFlipY(flipy); + flipper->SetFlipZ(flipz); + flipper->Update(); + peak_image = flipper->GetOutput(); + } + double rmse = 0.0; int iteration = 0; std::string name = "NOANCHOR"; // Load reference tractogram consisting of all known tracts std::vector< mitk::FiberBundle::Pointer > input_reference; mitk::FiberBundle::Pointer anchor_tractogram = mitk::IOUtil::Load(anchors_file); if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) ) { input_reference.push_back(anchor_tractogram); // Fit known tracts to peak image to obtain underexplained image MITK_INFO << "Fit anchor tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetTractograms(input_reference); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetPeakImage(peak_image); fitter->SetVerbose(true); fitter->SetMaskImage(mask); fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); rmse = fitter->GetRMSE(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[0] << " " << name << " RMSE: " << rmse << "\n"; name = ist::GetFilenameWithoutExtension(anchors_file); mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0); anchor_tracts->SetFiberColors(255,255,255); mitk::IOUtil::Save(anchor_tracts, out_folder + boost::lexical_cast((int)(100000*rms_diff[0])) + "_" + name + ".fib"); peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_" + name + ".nii.gz"); peak_image_writer->Update(); } if (use_weights || use_num_streamlines) { MITK_INFO << "Using tract weights as scores"; int c = 0; for (auto fib : input_candidates) { int mod = 1; double score = 0; if (use_weights) { score = fib->GetFiberWeight(0); mod = 100000; } else if (use_num_streamlines) score = fib->GetNumFibers(); fib->ColorFibersByOrientation(); std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(mod*score)) + "_" + bundle_name + ".fib"); float best_overlap = 0; int best_overlap_index = -1; int m_idx = 0; for (auto ref_mask : reference_masks) { float overlap = fib->GetOverlap(ref_mask); if (overlap>best_overlap) { best_overlap = overlap; best_overlap_index = m_idx; } ++m_idx; } unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } double weight_sum = 0; for (int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << "RMS_DIFF: " << setprecision(5) << score << " " << bundle_name << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; if (best_overlap_index>=0) logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(best_overlap_index)) << "\n"; else logfile << "No_overlap\n"; ++c; } } else if (!greedy_add) { MITK_INFO << "Fit candidate tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(true); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); fitter->SetTractograms(input_candidates); fitter->SetFitIndividualFibers(true); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); // vnl_vector log_rms_diff = rms_diff-rms_diff.min_value() + 1; // log_rms_diff = log_rms_diff.apply(std::log); // log_rms_diff /= log_rms_diff.max_value(); int c = 0; for (auto fib : input_candidates) { // fib->SetFiberWeights( log_rms_diff[c] ); // fib->ColorFibersByOrientation(); std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect if (filter_zero_weights) fib = fib->FilterByWeights(0); mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib"); float best_overlap = 0; int best_overlap_index = -1; int m_idx = 0; for (auto ref_mask : reference_masks) { float overlap = fib->GetOverlap(ref_mask); if (overlap>best_overlap) { best_overlap = overlap; best_overlap_index = m_idx; } ++m_idx; } unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } double weight_sum = 0; for (int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[c] << " " << bundle_name << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; if (best_overlap_index>=0) logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(best_overlap_index)) << "\n"; else logfile << "No_overlap\n"; ++c; } mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New(); out_fib = out_fib->AddBundles(input_candidates); out_fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib"); peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_AllCandidates.nii.gz"); peak_image_writer->Update(); } else { MITK_INFO << "RMSE: " << setprecision(5) << rmse; // fitter->SetPeakImage(peak_image); // Iteratively add candidate bundles in a greedy manner while (!input_candidates.empty()) { double next_rmse = rmse; double num_peaks = 0; mitk::FiberBundle::Pointer best_candidate = nullptr; PeakImgType::Pointer best_candidate_peak_image = nullptr; for (int i=0; i<(int)input_candidates.size(); ++i) { // WHY NECESSARY AGAIN?? itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(false); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); // ****************************** fitter->SetTractograms({input_candidates.at(i)}); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect fitter->Update(); std::cout.rdbuf (old); // <-- restore double candidate_rmse = fitter->GetRMSE(); if (candidate_rmseGetNumCoveredDirections(); best_candidate = fitter->GetTractograms().at(0); best_candidate_peak_image = fitter->GetUnderexplainedImage(); } } if (best_candidate.IsNull()) break; // fitter->SetPeakImage(peak_image); peak_image = best_candidate_peak_image; int i=0; std::vector< mitk::FiberBundle::Pointer > remaining_candidates; std::vector< std::string > remaining_candidate_files; for (auto fib : input_candidates) { if (fib!=best_candidate) { remaining_candidates.push_back(fib); remaining_candidate_files.push_back(candidate_tract_files.at(i)); } else name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i)); ++i; } input_candidates = remaining_candidates; candidate_tract_files = remaining_candidate_files; iteration++; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect // Save winning candidate if (filter_zero_weights) best_candidate = best_candidate->FilterByWeights(0); mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast(iteration) + "_" + name + ".fib"); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); peak_image_writer->Update(); // Calculate best overlap with reference masks for evaluation purposes float best_overlap = 0; int best_overlap_index = -1; i = 0; for (auto ref_mask : reference_masks) { float overlap = best_candidate->GetOverlap(ref_mask); if (overlap>best_overlap) { best_overlap = overlap; best_overlap_index = i; } ++i; } std::cout.rdbuf (old); // <-- restore logfile << "RMSE: " << setprecision(5) << rmse << " " << name << " " << num_peaks << "\n"; if (best_overlap_index>=0) logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(best_overlap_index)) << "\n"; else logfile << "No_overlap\n"; } } clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s"; 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/ExtractSimilarTracts.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp index 5cefcfc547..edf0a9c4a8 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp @@ -1,237 +1,237 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; 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()); } ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = mitk::IOUtil::Load(filename); 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:", "EU_MEAN (default), EU_STD, EU_MAX"); 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< 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); + mitk::FiberBundle::Pointer extracted = 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 { 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(); std::vector tracts = clusterer->GetOutCentroids(); ref_fib = mitk::FiberBundle::New(nullptr); ref_fib = ref_fib->AddBundles(tracts); mitk::IOUtil::Save(ref_fib, out_root + "centroids_" + ist::GetFilenameName(ref_bundle_files.at(c))); segmenter->SetInCentroids(ref_fib); } // 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))); + extracted = extracted->AddBundle(result); + mitk::IOUtil::Save(result, out_root + "extracted_" + 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 << "Extracted streamlines: " << extracted->GetNumFibers(); + mitk::IOUtil::Save(extracted, out_root + "extracted_streamlines.trk"); - MITK_INFO << "Streamlines remaining in candidate tractogram: " << fib->GetNumFibers(); - mitk::IOUtil::Save(fib, out_root + "candidate_tractogram.trk"); + MITK_INFO << "Residual streamlines: " << fib->GetNumFibers(); + mitk::IOUtil::Save(fib, out_root + "residual_streamlines.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/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitView.cpp index aa1703963f..1ca7905f97 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberFitView.cpp @@ -1,264 +1,265 @@ /*=================================================================== 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 "QmitkFiberFitView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkFiberFitView::VIEW_ID = "org.mitk.views.fiberfit"; using namespace mitk; QmitkFiberFitView::QmitkFiberFitView() : QmitkAbstractView() , m_Controls( nullptr ) { } // Destructor QmitkFiberFitView::~QmitkFiberFitView() { } void QmitkFiberFitView::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::QmitkFiberFitViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_StartButton, SIGNAL(clicked()), this, SLOT(StartFit()) ); connect( m_Controls->m_ImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DataSelectionChanged()) ); connect( m_Controls->m_TractBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DataSelectionChanged()) ); mitk::TNodePredicateDataType::Pointer isFib = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_TractBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TractBox->SetPredicate(isFib); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetPredicate( mitk::NodePredicateOr::New( mitk::NodePredicateAnd::New(isImage, is3D), mitk::TNodePredicateDataType::New()) ); DataSelectionChanged(); } } void QmitkFiberFitView::DataSelectionChanged() { if (m_Controls->m_TractBox->GetSelectedNode().IsNull() || m_Controls->m_ImageBox->GetSelectedNode().IsNull()) m_Controls->m_StartButton->setEnabled(false); else m_Controls->m_StartButton->setEnabled(true); } void QmitkFiberFitView::SetFocus() { DataSelectionChanged(); } void QmitkFiberFitView::StartFit() { if (m_Controls->m_TractBox->GetSelectedNode().IsNull() || m_Controls->m_ImageBox->GetSelectedNode().IsNull()) return; mitk::FiberBundle::Pointer input_tracts = dynamic_cast(m_Controls->m_TractBox->GetSelectedNode()->GetData()); mitk::DataNode::Pointer node = m_Controls->m_ImageBox->GetSelectedNode(); itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); mitk::Image::Pointer mitk_image = dynamic_cast(node->GetData()); mitk::PeakImage::Pointer mitk_peak_image = dynamic_cast(node->GetData()); if (mitk_peak_image.IsNotNull()) { typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitk_peak_image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer peak_image = caster->GetOutput(); fitter->SetPeakImage(peak_image); } else { if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { fitter->SetDiffImage(mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_image)); mitk::TensorModel<>* model = new mitk::TensorModel<>(); model->SetBvalue(1000); model->SetDiffusivity1(0.0010); model->SetDiffusivity2(0.00015); model->SetDiffusivity3(0.00015); model->SetGradientList(mitk::DiffusionPropertyHelper::GetGradientContainer(mitk_image)); fitter->SetSignalModel(model); } else { itk::FitFibersToImageFilter::DoubleImgType::Pointer scalar_image = itk::FitFibersToImageFilter::DoubleImgType::New(); mitk::CastToItkImage(mitk_image, scalar_image); fitter->SetScalarImage(scalar_image); } } if (m_Controls->m_ReguTypeBox->currentIndex()==0) fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (m_Controls->m_ReguTypeBox->currentIndex()==1) fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (m_Controls->m_ReguTypeBox->currentIndex()==2) fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (m_Controls->m_ReguTypeBox->currentIndex()==3) fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (m_Controls->m_ReguTypeBox->currentIndex()==4) fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->SetTractograms({input_tracts}); fitter->SetFitIndividualFibers(true); fitter->SetMaxIterations(20); fitter->SetVerbose(true); fitter->SetGradientTolerance(1e-5); fitter->SetLambda(m_Controls->m_ReguBox->value()); fitter->SetFilterOutliers(m_Controls->m_OutliersBox->isChecked()); fitter->Update(); mitk::FiberBundle::Pointer output_tracts = fitter->GetTractograms().at(0); + output_tracts->ColorFibersByFiberWeights(false, true); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(output_tracts); new_node->SetName("Fitted"); this->GetDataStorage()->Add(new_node, node); m_Controls->m_TractBox->GetSelectedNode()->SetVisibility(false); if (m_Controls->m_ResidualsBox->isChecked() && mitk_peak_image.IsNotNull()) { { mitk::PeakImage::ItkPeakImageType::Pointer itk_image = fitter->GetFittedImage(); mitk::Image::Pointer mitk_image = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(itk_image, mitk_image); mitk_image->SetVolume(itk_image->GetBufferPointer()); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(mitk_image); new_node->SetName("Fitted"); this->GetDataStorage()->Add(new_node, node); } { mitk::PeakImage::ItkPeakImageType::Pointer itk_image = fitter->GetResidualImage(); mitk::Image::Pointer mitk_image = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(itk_image, mitk_image); mitk_image->SetVolume(itk_image->GetBufferPointer()); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(mitk_image); new_node->SetName("Residual"); this->GetDataStorage()->Add(new_node, node); } { mitk::PeakImage::ItkPeakImageType::Pointer itk_image = fitter->GetUnderexplainedImage(); mitk::Image::Pointer mitk_image = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(itk_image, mitk_image); mitk_image->SetVolume(itk_image->GetBufferPointer()); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(mitk_image); new_node->SetName("Underexplained"); this->GetDataStorage()->Add(new_node, node); } { mitk::PeakImage::ItkPeakImageType::Pointer itk_image = fitter->GetOverexplainedImage(); mitk::Image::Pointer mitk_image = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(itk_image, mitk_image); mitk_image->SetVolume(itk_image->GetBufferPointer()); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(mitk_image); new_node->SetName("Overexplained"); this->GetDataStorage()->Add(new_node, node); } } else if (m_Controls->m_ResidualsBox->isChecked() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetFittedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(outImage); new_node->SetName("Fitted"); this->GetDataStorage()->Add(new_node, node); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetResidualImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(outImage); new_node->SetName("Residual"); this->GetDataStorage()->Add(new_node, node); } } else if (m_Controls->m_ResidualsBox->isChecked()) { { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetFittedImageScalar().GetPointer() ); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(outImage); new_node->SetName("Fitted"); this->GetDataStorage()->Add(new_node, node); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetResidualImageScalar().GetPointer() ); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData(outImage); new_node->SetName("Residual"); this->GetDataStorage()->Add(new_node, node); } } } void QmitkFiberFitView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { }