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 <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <itkImageDuplicator.h>
#include <boost/progress.hpp>
#define _USE_MATH_DEFINES
#include <math.h>
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; i<m_ImageSet->Size(); i++)
{
typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New();
duplicator->SetInputImage( m_ImageSet->GetElement(i) );
duplicator->Update();
set1->InsertElement(i, dynamic_cast<DirectionImageType*>(duplicator->GetOutput()));
}
for (unsigned int i=0; i<m_ReferenceImageSet->Size(); i++)
{
typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New();
duplicator->SetInputImage( m_ReferenceImageSet->GetElement(i) );
duplicator->Update();
set2->InsertElement(i, dynamic_cast<DirectionImageType*>(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<float>::max();
m_VarAngularError = 0.0;
m_AngularErrorVector.clear();
m_MeanLengthError = 0.0;
m_MedianLengthError = 0;
m_MaxLengthError = 0.0;
m_MinLengthError = itk::NumericTraits<float>::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<numReferenceImages; i++)
{
vnl_vector_fixed< PixelType, 3 > 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<numTestImages; i++)
{
vnl_vector_fixed< PixelType, 3 > 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<maxNumDirections; i++) // for each reference direction
{
bool missingDir = false;
vnl_vector_fixed< PixelType, 3 > refDir;
if (i<numRefDirs) // normalize if not null
refDir = refDirs.at(i);
else if (m_IgnoreMissingDirections)
continue;
else
missingDir = true;
for (int j=0; j<maxNumDirections; j++) // and each test direction
{
vnl_vector_fixed< PixelType, 3 > testDir;
if (j<numTestDirs) // normalize if not null
testDir = testDirs.at(j);
else if (m_IgnoreMissingDirections || missingDir)
continue;
else
missingDir = true;
// found missing direction
if (missingDir)
{
diffM[i][j] = -1;
continue;
}
// calculate angle between directions
diffM[i][j] = fabs(dot_product(refDir, testDir));
if (diffM[i][j] < maxAngularError)
maxAngularError = diffM[i][j];
if (diffM[i][j]>1.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<maxNumDirections; k++)
{
float error = -1;
int a,b; a=-1; b=-1;
bool missingDir = false;
// i: index of reference direction
// j: index of test direction
// find smalles error between two directions (large value -> small error)
for (int i=0; i<maxNumDirections; i++)
for (int j=0; j<maxNumDirections; j++)
{
if (diffM[i][j]>error && 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; i<numRefDirs; i++)
if (diffM_copy[i][b]>error)
{
error = diffM_copy[i][b];
a = i;
}
}
else if (b>=numTestDirs) // missing test direction (find next closest)
{
for (int i=0; i<numTestDirs; i++)
if (diffM_copy[a][i]>error)
{
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; i<m_AngularErrorVector.size(); i++)
{
float temp = m_AngularErrorVector.at(i) - m_MeanAngularError;
m_VarAngularError += temp*temp;
if ( m_AngularErrorVector.at(i)>m_MaxAngularError )
m_MaxAngularError = m_AngularErrorVector.at(i);
if ( m_AngularErrorVector.at(i)<m_MinAngularError )
m_MinAngularError = m_AngularErrorVector.at(i);
}
if (m_AngularErrorVector.size()>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; i<m_LengthErrorVector.size(); i++)
{
float temp = m_LengthErrorVector.at(i) - m_MeanLengthError;
m_VarLengthError += temp*temp;
if ( m_LengthErrorVector.at(i)>m_MaxLengthError )
m_MaxLengthError = m_LengthErrorVector.at(i);
if ( m_LengthErrorVector.at(i)<m_MinLengthError )
m_MinLengthError = m_LengthErrorVector.at(i);
}
if (m_LengthErrorVector.size()>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 <itkProcessObject.h>
#include <itkVectorContainer.h>
#include <itkImageSource.h>
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<Self> Pointer;
typedef SmartPointer<const Self> 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 <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <metaCommand.h>
#include "mitkCommandLineParser.h"
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <mitkIOUtil.h>
#include <iostream>
#include <fstream>
#include <itksys/SystemTools.hxx>
#include <mitkCoreObjectFactory.h>
#include <mitkFiberBundle.h>
#include <mitkTrackingForestHandler.h>
#define _USE_MATH_DEFINES
#include <math.h>
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<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
mitkCommandLineParser::StringContainerType imageFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["images"]);
- mitkCommandLineParser::StringContainerType wmMaskFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["wmmasks"]);
+ mitkCommandLineParser::StringContainerType wmMaskFiles;
+ if (parsedArgs.count("wmmasks"))
+ wmMaskFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["wmmasks"]);
mitkCommandLineParser::StringContainerType maskFiles;
if (parsedArgs.count("masks"))
maskFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["masks"]);
string forestFile = us::any_cast<string>(parsedArgs["forest"]);
mitkCommandLineParser::StringContainerType tractogramFiles;
if (parsedArgs.count("tractograms"))
tractogramFiles = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["tractograms"]);
int numTrees = 30;
if (parsedArgs.count("numtrees"))
numTrees = us::any_cast<int>(parsedArgs["numtrees"]);
int gmsamples = 50;
if (parsedArgs.count("gmsamples"))
gmsamples = us::any_cast<int>(parsedArgs["gmsamples"]);
float stepsize = -1;
if (parsedArgs.count("stepsize"))
stepsize = us::any_cast<float>(parsedArgs["stepsize"]);
int max_tree_depth = 50;
if (parsedArgs.count("max_tree_depth"))
max_tree_depth = us::any_cast<int>(parsedArgs["max_tree_depth"]);
double sample_fraction = 1.0;
if (parsedArgs.count("sample_fraction"))
sample_fraction = us::any_cast<float>(parsedArgs["sample_fraction"]);
MITK_INFO << "loading diffusion-weighted images";
std::vector< mitk::Image::Pointer > rawData;
for (unsigned int i=0; i<imageFiles.size(); i++)
{
mitk::Image::Pointer dwi = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadImage(imageFiles.at(i)).GetPointer());
rawData.push_back(dwi);
}
typedef itk::Image<unsigned char, 3> ItkUcharImgType;
MITK_INFO << "loading mask images";
std::vector< ItkUcharImgType::Pointer > maskImageVector;
for (unsigned int i=0; i<maskFiles.size(); i++)
{
mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(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<wmMaskFiles.size(); i++)
{
mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(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<tractogramFiles.size(); t++)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(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 <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <itkEvaluateDirectionImagesFilter.h>
#include <metaCommand.h>
#include "mitkCommandLineParser.h"
#include <itkTractsToVectorImageFilter.h>
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
#include <itksys/SystemTools.hxx>
#include <mitkCoreObjectFactory.h>
#define _USE_MATH_DEFINES
#include <math.h>
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<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
mitkCommandLineParser::StringContainerType referenceImages = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["reference"]);
mitkCommandLineParser::StringContainerType maskImages;
if (parsedArgs.count("mask"))
maskImages = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["mask"]);
string fibFile = us::any_cast<string>(parsedArgs["input"]);
float angularThreshold = 25;
if (parsedArgs.count("athresh"))
angularThreshold = us::any_cast<float>(parsedArgs["athresh"]);
float sizeThreshold = 0;
if (parsedArgs.count("sthresh"))
sizeThreshold = us::any_cast<float>(parsedArgs["sthresh"]);
int maxDirs = 0;
if (parsedArgs.count("maxdirs"))
maxDirs = us::any_cast<int>(parsedArgs["maxdirs"]);
string outRoot = us::any_cast<string>(parsedArgs["out"]);
bool verbose = false;
if (parsedArgs.count("verbose"))
verbose = us::any_cast<bool>(parsedArgs["verbose"]);
- bool ignore = false;
+ bool ignoreMissing = false;
if (parsedArgs.count("ignore"))
- ignore = us::any_cast<bool>(parsedArgs["ignore"]);
+ ignoreMissing = us::any_cast<bool>(parsedArgs["ignore"]);
+
+ bool ignoreEmpty = false;
+ if (parsedArgs.count("empty"))
+ ignoreEmpty = us::any_cast<bool>(parsedArgs["empty"]);
string fileID = "";
if (parsedArgs.count("fileID"))
fileID = us::any_cast<string>(parsedArgs["fileID"]);
try
{
typedef itk::Image<unsigned char, 3> 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::FiberBundle*>(mitk::IOUtil::LoadDataNode(fibFile)->GetData());
// load reference directions
ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New();
for (unsigned int i=0; i<referenceImages.size(); i++)
{
try
{
mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(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<float>::Pointer fOdfFilter = itk::TractsToVectorImageFilter<float>::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; i<directionImageContainer->Size(); i++)
{
itk::TractsToVectorImageFilter<float>::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i);
typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter<float>::ItkDirectionImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
string outfilename = outRoot;
outfilename.append("_DIRECTION_");
outfilename.append(boost::lexical_cast<string>(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<maskImages.size(); i++)
{
mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(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<string>(evaluationFilter->GetMeanAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMedianAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMaxAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMinAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(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<string>(evaluationFilter->GetMeanAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMedianAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMaxAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMinAngularError()));
sens.append(",");
sens.append(boost::lexical_cast<string>(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;
}