diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp index 56d48bbb46..a9d5174fe4 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp @@ -1,363 +1,372 @@ /*=================================================================== 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 __itkEvaluateDirectionImagesFilter_cpp #define __itkEvaluateDirectionImagesFilter_cpp #include "itkEvaluateDirectionImagesFilter.h" #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class PixelType > EvaluateDirectionImagesFilter< PixelType > ::EvaluateDirectionImagesFilter(): m_ImageSet(NULL), m_ReferenceImageSet(NULL), m_IgnoreMissingDirections(false), + m_IgnoreEmptyVoxels(false), m_Eps(0.0001) { this->SetNumberOfIndexedOutputs(2); } template< class PixelType > void EvaluateDirectionImagesFilter< PixelType >::GenerateData() { if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) return; DirectionImageContainerType::Pointer set1 = DirectionImageContainerType::New(); DirectionImageContainerType::Pointer set2 = DirectionImageContainerType::New(); for (unsigned int i=0; iSize(); i++) { typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); duplicator->SetInputImage( m_ImageSet->GetElement(i) ); duplicator->Update(); set1->InsertElement(i, dynamic_cast(duplicator->GetOutput())); } for (unsigned int i=0; iSize(); i++) { typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); duplicator->SetInputImage( m_ReferenceImageSet->GetElement(i) ); duplicator->Update(); set2->InsertElement(i, dynamic_cast(duplicator->GetOutput())); } m_ImageSet = set1; m_ReferenceImageSet = set2; // angular error image typename OutputImageType::Pointer outputImage = OutputImageType::New(); outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); this->SetNthOutput(0, outputImage); // length error image outputImage = OutputImageType::New(); outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); this->SetNthOutput(1, outputImage); if (m_MaskImage.IsNull()) { m_MaskImage = UCharImageType::New(); m_MaskImage->SetOrigin( outputImage->GetOrigin() ); m_MaskImage->SetRegions( outputImage->GetLargestPossibleRegion() ); m_MaskImage->SetSpacing( outputImage->GetSpacing() ); m_MaskImage->SetDirection( outputImage->GetDirection() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_MeanAngularError = 0.0; m_MedianAngularError = 0; m_MaxAngularError = 0.0; m_MinAngularError = itk::NumericTraits::max(); m_VarAngularError = 0.0; m_AngularErrorVector.clear(); m_MeanLengthError = 0.0; m_MedianLengthError = 0; m_MaxLengthError = 0.0; m_MinLengthError = itk::NumericTraits::max(); m_VarLengthError = 0.0; m_LengthErrorVector.clear(); if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) return; outputImage = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(0)); typename OutputImageType::Pointer outputImage2 = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); ImageRegionIterator< OutputImageType > oit2(outputImage2, outputImage2->GetLargestPossibleRegion()); ImageRegionIterator< UCharImageType > mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); int numTestImages = m_ImageSet->Size(); int numReferenceImages = m_ReferenceImageSet->Size(); int maxNumDirections = std::max(numReferenceImages, numTestImages); // matrix containing the angular error between the directions vnl_matrix< float > diffM; diffM.set_size(maxNumDirections, maxNumDirections); boost::progress_display disp(outputImage->GetLargestPossibleRegion().GetSize()[0]*outputImage->GetLargestPossibleRegion().GetSize()[1]*outputImage->GetLargestPossibleRegion().GetSize()[2]); while( !oit.IsAtEnd() ) { ++disp; if( mit.Get()!=1 ) { ++oit; ++oit2; ++mit; continue; } typename OutputImageType::IndexType index = oit.GetIndex(); float maxAngularError = 1.0; diffM.fill(10); // initialize with invalid error value // get number of valid directions (length > 0) int numRefDirs = 0; int numTestDirs = 0; std::vector< vnl_vector_fixed< PixelType, 3 > > testDirs; std::vector< vnl_vector_fixed< PixelType, 3 > > refDirs; for (int i=0; i refDir = m_ReferenceImageSet->GetElement(i)->GetPixel(index).GetVnlVector(); if (refDir.magnitude() > m_Eps ) { refDir.normalize(); refDirs.push_back(refDir); numRefDirs++; } } for (int i=0; i testDir = m_ImageSet->GetElement(i)->GetPixel(index).GetVnlVector(); if (testDir.magnitude() > m_Eps ) { testDir.normalize(); testDirs.push_back(testDir); numTestDirs++; } } + if (m_IgnoreEmptyVoxels && (numRefDirs==0 || numTestDirs==0) ) + { + ++oit; + ++oit2; + ++mit; + continue; + } + // i: index of reference direction // j: index of test direction for (int i=0; i refDir; if (i testDir; if (j1.0) diffM[i][j] = 1.0; } } float angularError = 0.0; float lengthError = 0.0; int counter = 0; vnl_matrix< float > diffM_copy = diffM; for (int k=0; k small error) for (int i=0; ierror && diffM[i][j]<2) // found valid error entry { error = diffM[i][j]; a = i; b = j; missingDir = false; } else if (diffM[i][j]<0 && error<0) // found missing direction { a = i; b = j; missingDir = true; } } if (a<0 || b<0 || (m_IgnoreMissingDirections && missingDir)) continue; // no more directions found if (a>=numRefDirs && b>=numTestDirs) { MITK_INFO << "ERROR: missing test and reference direction. should not be possible. check code."; continue; } // remove processed directions from error matrix diffM.set_row(a, 10.0); diffM.set_column(b, 10.0); if (a>=numRefDirs) // missing reference direction (find next closest) { for (int i=0; ierror) { error = diffM_copy[i][b]; a = i; } } else if (b>=numTestDirs) // missing test direction (find next closest) { for (int i=0; ierror) { error = diffM_copy[a][i]; b = i; } } float refLength = 0; float testLength = 1; if (a>=numRefDirs || b>=numTestDirs || error<0) error = 0; else { refLength = m_ReferenceImageSet->GetElement(a)->GetPixel(index).GetVnlVector().magnitude(); testLength = m_ImageSet->GetElement(b)->GetPixel(index).GetVnlVector().magnitude(); } m_LengthErrorVector.push_back( fabs(refLength-testLength) ); m_AngularErrorVector.push_back( acos(error)*180.0/M_PI ); m_MeanAngularError += m_AngularErrorVector.back(); m_MeanLengthError += m_LengthErrorVector.back(); angularError += m_AngularErrorVector.back(); lengthError += m_LengthErrorVector.back(); counter++; } if (counter>0) { lengthError /= counter; angularError /= counter; } oit2.Set(lengthError); oit.Set(angularError); ++oit; ++oit2; ++mit; } std::sort( m_AngularErrorVector.begin(), m_AngularErrorVector.end() ); m_MeanAngularError /= m_AngularErrorVector.size(); // mean for (unsigned int i=0; im_MaxAngularError ) m_MaxAngularError = m_AngularErrorVector.at(i); if ( m_AngularErrorVector.at(i)1) { m_VarAngularError /= (m_AngularErrorVector.size()-1); // variance // median if (m_AngularErrorVector.size()%2 == 0) m_MedianAngularError = 0.5*( m_AngularErrorVector.at( m_AngularErrorVector.size()/2 ) + m_AngularErrorVector.at( m_AngularErrorVector.size()/2+1 ) ); else m_MedianAngularError = m_AngularErrorVector.at( (m_AngularErrorVector.size()+1)/2 ) ; } std::sort( m_LengthErrorVector.begin(), m_LengthErrorVector.end() ); m_MeanLengthError /= m_LengthErrorVector.size(); // mean for (unsigned int i=0; im_MaxLengthError ) m_MaxLengthError = m_LengthErrorVector.at(i); if ( m_LengthErrorVector.at(i)1) { m_VarLengthError /= (m_LengthErrorVector.size()-1); // variance // median if (m_LengthErrorVector.size()%2 == 0) m_MedianLengthError = 0.5*( m_LengthErrorVector.at( m_LengthErrorVector.size()/2 ) + m_LengthErrorVector.at( m_LengthErrorVector.size()/2+1 ) ); else m_MedianLengthError = m_LengthErrorVector.at( (m_LengthErrorVector.size()+1)/2 ) ; } } } #endif // __itkEvaluateDirectionImagesFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h index 320d98de29..17a08617ca 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h @@ -1,112 +1,114 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkEvaluateDirectionImagesFilter_h_ #define __itkEvaluateDirectionImagesFilter_h_ #include #include #include namespace itk{ /** \brief Evaluates the voxel-wise angular error between two sets of directions. */ template< class PixelType > class EvaluateDirectionImagesFilter : public ImageSource< Image< PixelType, 3 > > { public: typedef EvaluateDirectionImagesFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageSource< Image< PixelType, 3 > > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef typename Superclass::OutputImageType OutputImageType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(EvaluateDirectionImagesFilter, ImageToImageFilter) typedef Vector< float, 3 > DirectionType; typedef Image< DirectionType, 3 > DirectionImageType; typedef VectorContainer< unsigned int, DirectionImageType::Pointer > DirectionImageContainerType; typedef Image< float, 3 > FloatImageType; typedef Image< bool, 3 > BoolImageType; typedef Image< unsigned char, 3 > UCharImageType; itkSetMacro( ImageSet , DirectionImageContainerType::Pointer) ///< test image containers itkSetMacro( ReferenceImageSet , DirectionImageContainerType::Pointer) ///< reference image containers itkSetMacro( MaskImage , UCharImageType::Pointer) ///< Calculation is only performed inside of the mask image. itkSetMacro( IgnoreMissingDirections , bool) ///< If in one voxel, the number of directions differs between the test container and the reference, the excess directions are ignored. Otherwise, the error to the next closest direction is calculated. + itkSetMacro( IgnoreEmptyVoxels , bool) ///< Don't increase error if either reference or test voxel is empty. /** Output statistics of the measured angular errors. */ itkGetMacro( MeanAngularError, float) itkGetMacro( MinAngularError, float) itkGetMacro( MaxAngularError, float) itkGetMacro( VarAngularError, float) itkGetMacro( MedianAngularError, float) /** Output statistics of the measured peak length errors. */ itkGetMacro( MeanLengthError, float) itkGetMacro( MinLengthError, float) itkGetMacro( MaxLengthError, float) itkGetMacro( VarLengthError, float) itkGetMacro( MedianLengthError, float) protected: EvaluateDirectionImagesFilter(); ~EvaluateDirectionImagesFilter() {} void GenerateData(); UCharImageType::Pointer m_MaskImage; DirectionImageContainerType::Pointer m_ImageSet; DirectionImageContainerType::Pointer m_ReferenceImageSet; bool m_IgnoreMissingDirections; + bool m_IgnoreEmptyVoxels; double m_MeanAngularError; double m_MedianAngularError; double m_MaxAngularError; double m_MinAngularError; double m_VarAngularError; std::vector< double > m_AngularErrorVector; double m_MeanLengthError; double m_MedianLengthError; double m_MaxLengthError; double m_MinLengthError; double m_VarLengthError; std::vector< double > m_LengthErrorVector; double m_Eps; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkEvaluateDirectionImagesFilter.cpp" #endif #endif //__itkEvaluateDirectionImagesFilter_h_ diff --git a/Modules/DiffusionImaging/MiniApps/DFTraining.cpp b/Modules/DiffusionImaging/MiniApps/DFTraining.cpp index a5b2e34103..62b0b18593 100755 --- a/Modules/DiffusionImaging/MiniApps/DFTraining.cpp +++ b/Modules/DiffusionImaging/MiniApps/DFTraining.cpp @@ -1,146 +1,148 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include int main(int argc, char* argv[]) { MITK_INFO << "DFTraining"; mitkCommandLineParser parser; parser.setTitle("Machine Learning Based Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("images", "i", mitkCommandLineParser::StringList, "DWIs:", "input diffusion-weighted images", us::Any(), false); - parser.addArgument("wmmasks", "w", mitkCommandLineParser::StringList, "WM-Masks:", "white matter mask images", us::Any(), false); + parser.addArgument("wmmasks", "w", mitkCommandLineParser::StringList, "WM-Masks:", "white matter mask images", us::Any()); parser.addArgument("tractograms", "t", mitkCommandLineParser::StringList, "Tractograms:", "input tractograms (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("masks", "m", mitkCommandLineParser::StringList, "Masks:", "mask images", us::Any()); parser.addArgument("forest", "f", mitkCommandLineParser::OutputFile, "Forest:", "output forest", us::Any(), false); parser.addArgument("stepsize", "s", mitkCommandLineParser::Float, "Stepsize:", "stepsize", us::Any()); parser.addArgument("gmsamples", "g", mitkCommandLineParser::Int, "Number of gray matter samples per voxel:", "Number of gray matter samples per voxel", us::Any()); parser.addArgument("numtrees", "n", mitkCommandLineParser::Int, "Number of trees:", "number of trees", us::Any()); parser.addArgument("max_tree_depth", "d", mitkCommandLineParser::Int, "Max. tree depth:", "maximum tree depth", us::Any()); parser.addArgument("sample_fraction", "sf", mitkCommandLineParser::Float, "Sample fraction:", "fraction of samples used per tree", us::Any()); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType imageFiles = us::any_cast(parsedArgs["images"]); - mitkCommandLineParser::StringContainerType wmMaskFiles = us::any_cast(parsedArgs["wmmasks"]); + mitkCommandLineParser::StringContainerType wmMaskFiles; + if (parsedArgs.count("wmmasks")) + wmMaskFiles = us::any_cast(parsedArgs["wmmasks"]); mitkCommandLineParser::StringContainerType maskFiles; if (parsedArgs.count("masks")) maskFiles = us::any_cast(parsedArgs["masks"]); string forestFile = us::any_cast(parsedArgs["forest"]); mitkCommandLineParser::StringContainerType tractogramFiles; if (parsedArgs.count("tractograms")) tractogramFiles = us::any_cast(parsedArgs["tractograms"]); int numTrees = 30; if (parsedArgs.count("numtrees")) numTrees = us::any_cast(parsedArgs["numtrees"]); int gmsamples = 50; if (parsedArgs.count("gmsamples")) gmsamples = us::any_cast(parsedArgs["gmsamples"]); float stepsize = -1; if (parsedArgs.count("stepsize")) stepsize = us::any_cast(parsedArgs["stepsize"]); int max_tree_depth = 50; if (parsedArgs.count("max_tree_depth")) max_tree_depth = us::any_cast(parsedArgs["max_tree_depth"]); double sample_fraction = 1.0; if (parsedArgs.count("sample_fraction")) sample_fraction = us::any_cast(parsedArgs["sample_fraction"]); MITK_INFO << "loading diffusion-weighted images"; std::vector< mitk::Image::Pointer > rawData; for (unsigned int i=0; i(mitk::IOUtil::LoadImage(imageFiles.at(i)).GetPointer()); rawData.push_back(dwi); } typedef itk::Image ItkUcharImgType; MITK_INFO << "loading mask images"; std::vector< ItkUcharImgType::Pointer > maskImageVector; for (unsigned int i=0; i(mitk::IOUtil::LoadImage(maskFiles.at(i)).GetPointer()); ItkUcharImgType::Pointer mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); maskImageVector.push_back(mask); } MITK_INFO << "loading white matter mask images"; std::vector< ItkUcharImgType::Pointer > wmMaskImageVector; for (unsigned int i=0; i(mitk::IOUtil::LoadImage(wmMaskFiles.at(i)).GetPointer()); ItkUcharImgType::Pointer wmmask = ItkUcharImgType::New(); mitk::CastToItkImage(img, wmmask); wmMaskImageVector.push_back(wmmask); } MITK_INFO << "loading tractograms"; std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int t=0; t(mitk::IOUtil::Load(tractogramFiles.at(t)).at(0).GetPointer()); tractograms.push_back(fib); } mitk::TrackingForestHandler<> forestHandler; forestHandler.SetRawData(rawData); forestHandler.SetTractograms(tractograms); forestHandler.SetNumTrees(numTrees); forestHandler.SetMaxTreeDepth(max_tree_depth); forestHandler.SetGrayMatterSamplesPerVoxel(gmsamples); forestHandler.SetSampleFraction(sample_fraction); forestHandler.SetStepSize(stepsize); forestHandler.StartTraining(); forestHandler.SaveForest(forestFile); return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp index 2154df23e4..64d7707814 100755 --- a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp +++ b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp @@ -1,312 +1,319 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Local Directional Fiber Plausibility"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(" "); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("reference", "r", mitkCommandLineParser::StringList, "Reference images:", "reference direction images", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::StringList, "Masks:", "mask images"); parser.addArgument("athresh", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); parser.addArgument("sthresh", "s", mitkCommandLineParser::Float, "Size threshold:", "Relative peak size threshold per voxel.", 0.0, true); parser.addArgument("maxdirs", "md", mitkCommandLineParser::Int, "Max. Clusters:", "Maximum number of fiber clusters.", 0, true); parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results"); parser.addArgument("ignore", "n", mitkCommandLineParser::Bool, "Ignore:", "don't increase error for missing or too many directions"); + parser.addArgument("empty", "e", mitkCommandLineParser::Bool, "Empty Voxels:", "don't increase error for empty voxels"); parser.addArgument("fileID", "id", mitkCommandLineParser::String, "ID:", "optional ID field"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType referenceImages = us::any_cast(parsedArgs["reference"]); mitkCommandLineParser::StringContainerType maskImages; if (parsedArgs.count("mask")) maskImages = us::any_cast(parsedArgs["mask"]); string fibFile = us::any_cast(parsedArgs["input"]); float angularThreshold = 25; if (parsedArgs.count("athresh")) angularThreshold = us::any_cast(parsedArgs["athresh"]); float sizeThreshold = 0; if (parsedArgs.count("sthresh")) sizeThreshold = us::any_cast(parsedArgs["sthresh"]); int maxDirs = 0; if (parsedArgs.count("maxdirs")) maxDirs = us::any_cast(parsedArgs["maxdirs"]); string outRoot = us::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); - bool ignore = false; + bool ignoreMissing = false; if (parsedArgs.count("ignore")) - ignore = us::any_cast(parsedArgs["ignore"]); + ignoreMissing = us::any_cast(parsedArgs["ignore"]); + + bool ignoreEmpty = false; + if (parsedArgs.count("empty")) + ignoreEmpty = us::any_cast(parsedArgs["empty"]); string fileID = ""; if (parsedArgs.count("fileID")) fileID = us::any_cast(parsedArgs["fileID"]); try { typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); // load reference directions ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); } catch(...){ std::cout << "could not load: " << referenceImages.at(i); } } ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); itkMaskImage->SetSpacing( dirImg->GetSpacing() ); itkMaskImage->SetOrigin( dirImg->GetOrigin() ); itkMaskImage->SetDirection( dirImg->GetDirection() ); itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->Allocate(); itkMaskImage->FillBuffer(1); // extract directions from fiber bundle itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); fOdfFilter->SetNormalizeVectors(true); fOdfFilter->SetUseWorkingCopy(false); fOdfFilter->SetSizeThreshold(sizeThreshold); fOdfFilter->SetMaxNumDirections(maxDirs); fOdfFilter->Update(); ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); if (verbose) { // write vector field mitk::FiberBundle::Pointer directions = fOdfFilter->GetOutputFiberBundle(); string outfilename = outRoot; outfilename.append("_VECTOR_FIELD.fib"); mitk::IOUtil::SaveBaseData(directions.GetPointer(), outfilename ); // write direction images for (unsigned int i=0; iSize(); i++) { itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_DIRECTION_"); outfilename.append(boost::lexical_cast(i)); outfilename.append(".nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(itkImg); writer->Update(); } // write num direction image { ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_NUM_DIRECTIONS.nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(numDirImage); writer->Update(); } } string logFile = outRoot; logFile.append("_ANGULAR_ERROR.csv"); ofstream file; file.open (logFile.c_str()); if (maskImages.size()>0) { for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(maskImages.at(i))->GetData()); mitk::CastToItkImage(mitkMaskImage, itkMaskImage); // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); - evaluationFilter->SetIgnoreMissingDirections(ignore); + evaluationFilter->SetIgnoreMissingDirections(ignoreMissing); + evaluationFilter->SetIgnoreEmptyVoxels(ignoreEmpty); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string maskFileName = itksys::SystemTools::GetFilenameWithoutExtension(maskImages.at(i)); unsigned found = maskFileName.find_last_of("_"); string sens = itksys::SystemTools::GetFilenameWithoutLastExtension(fibFile); if (!fileID.empty()) sens = fileID; sens.append(","); sens.append(maskFileName.substr(found+1)); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; } } else { // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); - evaluationFilter->SetIgnoreMissingDirections(ignore); + evaluationFilter->SetIgnoreMissingDirections(ignoreMissing); + evaluationFilter->SetIgnoreEmptyVoxels(ignoreEmpty); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string sens = itksys::SystemTools::GetFilenameWithoutLastExtension(fibFile); if (!fileID.empty()) sens = fileID; sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; } file.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; }