diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt
index 5f79886390..6a26aa8477 100755
--- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt
+++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt
@@ -1,64 +1,65 @@
option(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF)
if(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS)
# needed include directories
include_directories(
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_BINARY_DIR}
)
project( MitkDiffusionMiniApps )
# fill in the standalone executables here
set(DIFFUSIONMINIAPPS
mitkDiffusionMiniApps
)
# set additional files here
set(DIFFUSIONCORE_ADDITIONAL_FILES
MiniAppManager.cpp
FileFormatConverter.cpp
TensorReconstruction.cpp
QballReconstruction.cpp
DiffusionIndices.cpp
CopyGeometry.cpp
GibbsTracking.cpp
StreamlineTracking.cpp
FiberProcessing.cpp
LocalDirectionalFiberPlausibility.cpp
#TractogramAngularError.cpp
FiberDirectionExtraction.cpp
PeakExtraction.cpp
PeaksAngularError.cpp
MultishellMethods.cpp
Fiberfox.cpp
ExportShImage.cpp
+ ImportShImage.cpp
NetworkCreation.cpp
NetworkStatistics.cpp
DwiDenoising.cpp
)
# deprecated
# FOREACH(tool ${DIFFUSIONMINIAPPS})
# ADD_EXECUTABLE(
# ${tool}
# ${tool}.cpp
# ${DIFFUSIONCORE_ADDITIONAL_FILES}
# )
# TARGET_LINK_LIBRARIES(
# ${tool}
# ${ALL_LIBRARIES} )
# ENDFOREACH(tool)
mitk_create_executable(DiffusionMiniApps
DEPENDS MitkDiffusionCore MitkFiberTracking MitkConnectomics
PACKAGE_DEPENDS ITK|ITKDiffusionTensorImage
)
if(EXECUTABLE_IS_ENABLED)
MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET})
endif()
endif()
diff --git a/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp b/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp
index 5f9c51eba6..5495b3093e 100755
--- a/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp
+++ b/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp
@@ -1,97 +1,93 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <mitkDiffusionImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkIOUtil.h>
#include <mitkNrrdDiffusionImageWriter.h>
#include "ctkCommandLineParser.h"
using namespace mitk;
#include "ctkCommandLineParser.h"
int CopyGeometry(int argc, char* argv[])
{
ctkCommandLineParser parser;
parser.setTitle("Copy Geometry");
parser.setCategory("Preprocessing Tools");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("in", "i", ctkCommandLineParser::InputFile, "Input:", "input image", us::Any(), false);
parser.addArgument("ref", "r", ctkCommandLineParser::InputFile, "Reference:", "reference image", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::OutputFile, "Output:", "output image", us::Any(), false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
// mandatory arguments
string imageName = us::any_cast<string>(parsedArgs["in"]);
string refImage = us::any_cast<string>(parsedArgs["ref"]);
string outImage = us::any_cast<string>(parsedArgs["out"]);
try
{
-
- MITK_INFO << "Loading image " << imageName;
const std::string s1="", s2="";
std::vector<BaseData::Pointer> infile = BaseDataIO::LoadBaseDataFromFile( refImage, s1, s2, false );
Image::Pointer source = dynamic_cast<Image*>(infile.at(0).GetPointer());
infile = BaseDataIO::LoadBaseDataFromFile( imageName, s1, s2, false );
Image::Pointer target = dynamic_cast<Image*>(infile.at(0).GetPointer());
mitk::BaseGeometry* s_geom = source->GetGeometry();
mitk::BaseGeometry* t_geom = target->GetGeometry();
t_geom->SetIndexToWorldTransform(s_geom->GetIndexToWorldTransform());
target->SetGeometry(t_geom);
if ( dynamic_cast<DiffusionImage<short>*>(target.GetPointer()) )
{
- MITK_INFO << "Writing " << outImage;
DiffusionImage<short>::Pointer dwi = dynamic_cast<DiffusionImage<short>*>(target.GetPointer());
NrrdDiffusionImageWriter<short>::Pointer writer = NrrdDiffusionImageWriter<short>::New();
writer->SetFileName(outImage);
writer->SetInput(dwi);
writer->Update();
}
else
mitk::IOUtil::SaveImage(target, outImage);
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
- MITK_INFO << "DONE";
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(CopyGeometry);
diff --git a/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp b/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp
index 292e94f8f3..8e768840f7 100644
--- a/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp
+++ b/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp
@@ -1,152 +1,149 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <itkExceptionObject.h>
#include <itkImageFileWriter.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkQBallImage.h>
#include <itkTensorDerivedMeasurementsFilter.h>
#include <itkDiffusionQballGeneralizedFaImageFilter.h>
#include <mitkTensorImage.h>
#include "ctkCommandLineParser.h"
#include <boost/algorithm/string.hpp>
#include <itksys/SystemTools.hxx>
#include <itkMultiThreader.h>
/**
* Calculate indices derived from Qball or tensor images
*/
int DiffusionIndices(int argc, char* argv[])
{
+ MITK_INFO << "DiffusionIndices";
ctkCommandLineParser parser;
parser.setTitle("Diffusion Indices");
parser.setCategory("Diffusion Related Measures");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input:", "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", us::Any(), false);
parser.addArgument("index", "idx", ctkCommandLineParser::String, "Index:", "index (fa, gfa, ra, ad, rd, ca, l2, l3, md)", us::Any(), false);
parser.addArgument("outFile", "o", ctkCommandLineParser::OutputFile, "Output:", "output file", us::Any(), false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string inFileName = us::any_cast<string>(parsedArgs["input"]);
string index = us::any_cast<string>(parsedArgs["index"]);
string outFileName = us::any_cast<string>(parsedArgs["outFile"]);
string ext = itksys::SystemTools::GetFilenameLastExtension(outFileName);
if (ext.empty())
outFileName += ".nrrd";
try
{
// load input image
const std::string s1="", s2="";
std::vector<mitk::BaseData::Pointer> infile = mitk::BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false );
if( boost::algorithm::ends_with(inFileName, ".qbi") && index=="gfa" )
{
typedef itk::Vector<float, QBALL_ODFSIZE> OdfVectorType;
typedef itk::Image<OdfVectorType,3> ItkQballImageType;
mitk::QBallImage::Pointer mitkQballImage = dynamic_cast<mitk::QBallImage*>(infile.at(0).GetPointer());
ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New();
mitk::CastToItkImage(mitkQballImage, itk_qbi);
typedef itk::DiffusionQballGeneralizedFaImageFilter<float,float,QBALL_ODFSIZE> GfaFilterType;
GfaFilterType::Pointer gfaFilter = GfaFilterType::New();
gfaFilter->SetInput(itk_qbi);
gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD);
gfaFilter->Update();
- MITK_INFO << "Writing " << outFileName;
-
itk::ImageFileWriter< itk::Image<float,3> >::Pointer fileWriter = itk::ImageFileWriter< itk::Image<float,3> >::New();
fileWriter->SetInput(gfaFilter->GetOutput());
fileWriter->SetFileName(outFileName);
fileWriter->Update();
}
else if( boost::algorithm::ends_with(inFileName, ".dti") )
{
typedef itk::Image< itk::DiffusionTensor3D<float>, 3 > ItkTensorImage;
mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast<mitk::TensorImage*>(infile.at(0).GetPointer());
ItkTensorImage::Pointer itk_dti = ItkTensorImage::New();
mitk::CastToItkImage(mitkTensorImage, itk_dti);
typedef itk::TensorDerivedMeasurementsFilter<float> MeasurementsType;
MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New();
measurementsCalculator->SetInput(itk_dti.GetPointer() );
if(index=="fa")
measurementsCalculator->SetMeasure(MeasurementsType::FA);
else if(index=="ra")
measurementsCalculator->SetMeasure(MeasurementsType::RA);
else if(index=="ad")
measurementsCalculator->SetMeasure(MeasurementsType::AD);
else if(index=="rd")
measurementsCalculator->SetMeasure(MeasurementsType::RD);
else if(index=="ca")
measurementsCalculator->SetMeasure(MeasurementsType::CA);
else if(index=="l2")
measurementsCalculator->SetMeasure(MeasurementsType::L2);
else if(index=="l3")
measurementsCalculator->SetMeasure(MeasurementsType::L3);
else if(index=="md")
measurementsCalculator->SetMeasure(MeasurementsType::MD);
else
{
MITK_WARN << "No valid diffusion index for input image (tensor image) defined";
return EXIT_FAILURE;
}
measurementsCalculator->Update();
- MITK_INFO << "Writing " << outFileName;
-
itk::ImageFileWriter< itk::Image<float,3> >::Pointer fileWriter = itk::ImageFileWriter< itk::Image<float,3> >::New();
fileWriter->SetInput(measurementsCalculator->GetOutput());
fileWriter->SetFileName(outFileName);
fileWriter->Update();
}
else
MITK_INFO << "Diffusion index " << index << " not supported for supplied file type.";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(DiffusionIndices);
diff --git a/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp b/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp
index 6e7ecdc5d0..28f0d0402d 100644
--- a/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp
+++ b/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp
@@ -1,168 +1,165 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <mitkBaseDataIOFactory.h>
#include "ctkCommandLineParser.h"
#include <boost/algorithm/string.hpp>
#include <DiffusionWeightedImages/mitkDiffusionImage.h>
#include <itkNonLocalMeansDenoisingFilter.h>
#include <mitkNrrdDiffusionImageWriter.h>
#include <itkImage.h>
typedef mitk::DiffusionImage<short> DiffusionImageType;
typedef itk::Image<short, 3> ImageType;
mitk::BaseData::Pointer LoadFile(std::string filename)
{
if( filename.empty() )
return NULL;
const std::string s1="", s2="";
std::vector<mitk::BaseData::Pointer> infile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false );
if( infile.empty() )
{
MITK_INFO << "File " << filename << " could not be read!";
return NULL;
}
mitk::BaseData::Pointer baseData = infile.at(0);
return baseData;
}
/**
* Denoises DWI using the Nonlocal - Means algorithm
*/
int DwiDenoising(int argc, char* argv[])
{
+ MITK_INFO << "DwiDenoising";
ctkCommandLineParser parser;
parser.setTitle("DWI Denoising");
parser.setCategory("Preprocessing Tools");
parser.setContributor("MBI");
parser.setDescription("Denoising for diffusion weighted images using a non-local means algorithm.");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input:", "input image (DWI)", us::Any(), false);
parser.addArgument("variance", "v", ctkCommandLineParser::Float, "Variance:", "noise variance", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::InputFile, "Mask:", "brainmask for input image", us::Any(), true);
parser.addArgument("search", "s", ctkCommandLineParser::Int, "Search radius:", "search radius", us::Any(), true);
parser.addArgument("compare", "c", ctkCommandLineParser::Int, "Comparison radius:", "comparison radius", us::Any(), true);
parser.addArgument("joint", "j", ctkCommandLineParser::Bool, "Joint information:", "use joint information");
parser.addArgument("rician", "r", ctkCommandLineParser::Bool, "Rician adaption:", "use rician adaption");
parser.changeParameterGroup("Output", "Output of this miniapp");
parser.addArgument("output", "o", ctkCommandLineParser::OutputFile, "Output:", "output image (DWI)", us::Any(), false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string inFileName = us::any_cast<string>(parsedArgs["input"]);
double variance = static_cast<double>(us::any_cast<float>(parsedArgs["variance"]));
string maskName;
if (parsedArgs.count("mask"))
maskName = us::any_cast<string>(parsedArgs["mask"]);
string outFileName = us::any_cast<string>(parsedArgs["output"]);
// boost::algorithm::erase_all(outFileName, ".dwi");
int search = 4;
if (parsedArgs.count("search"))
search = us::any_cast<int>(parsedArgs["search"]);
int compare = 1;
if (parsedArgs.count("compare"))
compare = us::any_cast<int>(parsedArgs["compare"]);
bool joint = false;
if (parsedArgs.count("joint"))
joint = true;
bool rician = false;
if (parsedArgs.count("rician"))
rician = true;
try
{
if( boost::algorithm::ends_with(inFileName, ".dwi"))
{
DiffusionImageType::Pointer dwi = dynamic_cast<DiffusionImageType*>(LoadFile(inFileName).GetPointer());
itk::NonLocalMeansDenoisingFilter<short>::Pointer filter = itk::NonLocalMeansDenoisingFilter<short>::New();
filter->SetNumberOfThreads(12);
filter->SetInputImage(dwi->GetVectorImage());
if (!maskName.empty())
{
mitk::Image::Pointer mask = dynamic_cast<mitk::Image*>(LoadFile(maskName).GetPointer());
ImageType::Pointer itkMask = ImageType::New();
mitk::CastToItkImage(mask, itkMask);
filter->SetInputMask(itkMask);
}
filter->SetUseJointInformation(joint);
filter->SetUseRicianAdaption(rician);
filter->SetSearchRadius(search);
filter->SetComparisonRadius(compare);
filter->SetVariance(variance);
filter->Update();
DiffusionImageType::Pointer output = DiffusionImageType::New();
output->SetVectorImage(filter->GetOutput());
output->SetReferenceBValue(dwi->GetReferenceBValue());
output->SetDirections(dwi->GetDirections());
output->InitializeFromVectorImage();
// std::stringstream name;
// name << outFileName << "_NLM_" << search << "-" << compare << "-" << variance << ".dwi";
- MITK_INFO << "Writing: " << outFileName;
-
mitk::NrrdDiffusionImageWriter<short>::Pointer writer = mitk::NrrdDiffusionImageWriter<short>::New();
writer->SetInput(output);
writer->SetFileName(outFileName/*.str()*/);
writer->Update();
-
- MITK_INFO << "Finish!";
}
else
{
MITK_INFO << "Only supported for .dwi!";
}
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(DwiDenoising);
diff --git a/Modules/DiffusionImaging/MiniApps/ExportShImage.cpp b/Modules/DiffusionImaging/MiniApps/ExportShImage.cpp
index e81690ffb5..8bbd36c37b 100755
--- a/Modules/DiffusionImaging/MiniApps/ExportShImage.cpp
+++ b/Modules/DiffusionImaging/MiniApps/ExportShImage.cpp
@@ -1,140 +1,138 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <metaCommand.h>
#include "ctkCommandLineParser.h"
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <itkImageFileReader.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <itkShCoefficientImageExporter.h>
#include <itkFlipImageFilter.h>
#define _USE_MATH_DEFINES
#include <math.h>
template<int shOrder>
int StartShConversion(int argc, char* argv[])
{
ctkCommandLineParser parser;
parser.setTitle("Export SH Image");
parser.setCategory("Preprocessing Tools");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input:", "MITK SH image", us::Any(), false);
parser.addArgument("output", "o", ctkCommandLineParser::InputFile, "Output", "MRtrix SH image", us::Any(), false);
parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "SH order:", "spherical harmonics order");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string inFile = us::any_cast<string>(parsedArgs["input"]);
string outFile = us::any_cast<string>(parsedArgs["output"]);
try
{
typedef itk::Image< float, 4 > OutImageType;
typedef itk::Image< itk::Vector< float, (shOrder*shOrder + shOrder + 2)/2 + shOrder >, 3 > InputImageType;
typename InputImageType::Pointer itkInImage = InputImageType::New();
typedef itk::ImageFileReader< InputImageType > ReaderType;
typename ReaderType::Pointer reader = ReaderType::New();
MITK_INFO << "reading " << inFile;
reader->SetFileName(inFile.c_str());
reader->Update();
itkInImage = reader->GetOutput();
// extract directions from fiber bundle
typename itk::ShCoefficientImageExporter<float, shOrder>::Pointer filter = itk::ShCoefficientImageExporter<float,shOrder>::New();
filter->SetInputImage(itkInImage);
filter->GenerateData();
OutImageType::Pointer outImage = filter->GetOutputImage();
typedef itk::ImageFileWriter< OutImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
- MITK_INFO << "writing " << outFile;
writer->SetFileName(outFile.c_str());
writer->SetInput(outImage);
writer->Update();
-
- MITK_INFO << "DONE";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
int ExportShImage(int argc, char* argv[])
{
+ MITK_INFO << "ExportShImage";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input image", "MITK SH image", us::Any(), false);
parser.addArgument("output", "o", ctkCommandLineParser::OutputFile, "Output image", "MRtrix SH image", us::Any(), false);
parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order");
parser.setCategory("Preprocessing Tools");
parser.setTitle("Export SH Image");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
int shOrder = -1;
if (parsedArgs.count("shOrder"))
shOrder = us::any_cast<int>(parsedArgs["shOrder"]);
switch (shOrder)
{
case 4:
return StartShConversion<4>(argc, argv);
case 6:
return StartShConversion<6>(argc, argv);
case 8:
return StartShConversion<8>(argc, argv);
case 10:
return StartShConversion<10>(argc, argv);
case 12:
return StartShConversion<12>(argc, argv);
}
return EXIT_FAILURE;
}
RegisterDiffusionMiniApp(ExportShImage);
diff --git a/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp b/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
index e6cb2efd48..1488960489 100755
--- a/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
@@ -1,185 +1,182 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <metaCommand.h>
#include "ctkCommandLineParser.h"
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <itkEvaluateDirectionImagesFilter.h>
#include <itkTractsToVectorImageFilter.h>
#include <mitkCoreObjectFactory.h>
#define _USE_MATH_DEFINES
#include <math.h>
int FiberDirectionExtraction(int argc, char* argv[])
{
+ MITK_INFO << "FiberDirectionExtraction";
ctkCommandLineParser parser;
parser.setTitle("Fiber Direction Extraction");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::InputFile, "Mask:", "mask image");
parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true);
parser.addArgument("peakthresh", "t", ctkCommandLineParser::Float, "Peak size threshold:", "peak size threshold relative to largest peak in voxel", 0.2, true);
parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results");
parser.addArgument("numdirs", "d", ctkCommandLineParser::Int, "Max. num. directions:", "maximum number of fibers per voxel", 3, true);
parser.addArgument("normalize", "n", ctkCommandLineParser::Bool, "Normalize:", "normalize vectors");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string fibFile = us::any_cast<string>(parsedArgs["input"]);
string maskImage("");
if (parsedArgs.count("mask"))
maskImage = us::any_cast<string>(parsedArgs["mask"]);
float peakThreshold = 0.2;
if (parsedArgs.count("peakthresh"))
peakThreshold = us::any_cast<float>(parsedArgs["peakthresh"]);
float angularThreshold = 25;
if (parsedArgs.count("athresh"))
angularThreshold = us::any_cast<float>(parsedArgs["athresh"]);
string outRoot = us::any_cast<string>(parsedArgs["out"]);
bool verbose = false;
if (parsedArgs.count("verbose"))
verbose = us::any_cast<bool>(parsedArgs["verbose"]);
int maxNumDirs = 3;
if (parsedArgs.count("numdirs"))
maxNumDirs = us::any_cast<int>(parsedArgs["numdirs"]);
bool normalize = false;
if (parsedArgs.count("normalize"))
normalize = us::any_cast<bool>(parsedArgs["normalize"]);
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;
// load fiber bundle
mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(mitk::IOUtil::LoadDataNode(fibFile)->GetData());
// load/create mask image
ItkUcharImgType::Pointer itkMaskImage = NULL;
if (maskImage.compare("")!=0)
{
MITK_INFO << "Using mask image";
itkMaskImage = ItkUcharImgType::New();
mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(maskImage)->GetData());
mitk::CastToItkImage<ItkUcharImgType>(mitkMaskImage, itkMaskImage);
}
// 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(normalize);
fOdfFilter->SetUseWorkingCopy(false);
fOdfFilter->SetSizeThreshold(peakThreshold);
fOdfFilter->SetMaxNumDirections(maxNumDirs);
fOdfFilter->Update();
ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer();
// 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");
- MITK_INFO << "writing " << outfilename;
writer->SetFileName(outfilename.c_str());
writer->SetInput(itkImg);
writer->Update();
}
if (verbose)
{
// write vector field
mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle();
mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters();
for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it)
{
if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) {
string outfilename = outRoot;
outfilename.append("_VECTOR_FIELD.fib");
(*it)->SetFileName( outfilename.c_str() );
(*it)->DoWrite( directions.GetPointer() );
}
}
// 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");
- MITK_INFO << "writing " << outfilename;
writer->SetFileName(outfilename.c_str());
writer->SetInput(numDirImage);
writer->Update();
}
}
-
- MITK_INFO << "DONE";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(FiberDirectionExtraction);
diff --git a/Modules/DiffusionImaging/MiniApps/FiberExtraction.cpp b/Modules/DiffusionImaging/MiniApps/FiberExtraction.cpp
index 6b536d7874..cab8803780 100755
--- a/Modules/DiffusionImaging/MiniApps/FiberExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/FiberExtraction.cpp
@@ -1,159 +1,158 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkBaseDataIOFactory.h>
#include <metaCommand.h>
#include "ctkCommandLineParser.h"
#include <usAny.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <mitkCoreObjectFactory.h>
#include <mitkPlanarFigure.h>
#include <mitkPlanarFigureComposite.h>
#include <mitkFiberBundleX.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#define _USE_MATH_DEFINES
#include <math.h>
int FiberExtraction(int argc, char* argv[])
{
+ MITK_INFO << "FiberExtraction";
ctkCommandLineParser parser;
parser.setTitle("Fiber Extraction");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setContributor("MBI");
parser.setDescription("");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::String, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false);
parser.addArgument("planfirgure1", "pf1", ctkCommandLineParser::String, "Figure 1:", "first ROI", us::Any(), false);
parser.addArgument("planfirgure2", "pf2", ctkCommandLineParser::String, "Figure 2:", "second ROI", us::Any());
parser.addArgument("operation", "op", ctkCommandLineParser::String, "Operation:", "logical operation (AND, OR, NOT)", us::Any());
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string inFib = us::any_cast<string>(parsedArgs["input"]);
string outFib = us::any_cast<string>(parsedArgs["out"]);
string pf1_path = us::any_cast<string>(parsedArgs["planfirgure1"]);
string operation("");
string pf2_path("");
if (parsedArgs.count("operation"))
{
operation = us::any_cast<string>(parsedArgs["operation"]);
if (parsedArgs.count("planfirgure2") && (operation=="AND" || operation=="OR"))
pf2_path = us::any_cast<string>(parsedArgs["planfirgure2"]);
}
try
{
typedef itk::Image<unsigned char, 3> ItkUcharImgType;
// load fiber bundle
mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(mitk::IOUtil::LoadDataNode(inFib)->GetData());
mitk::FiberBundleX::Pointer result;
mitk::BaseData::Pointer input1 = mitk::IOUtil::LoadDataNode(pf1_path)->GetData();
mitk::PlanarFigure::Pointer pf1 = dynamic_cast<mitk::PlanarFigure*>(input1.GetPointer());
if (pf1.IsNotNull())
{
mitk::BaseData::Pointer input2;
mitk::PlanarFigure::Pointer pf2;
if (!pf2_path.empty())
{
input2 = mitk::IOUtil::LoadDataNode(pf2_path)->GetData();
pf2 = dynamic_cast<mitk::PlanarFigure*>(input2.GetPointer());
}
mitk::PlanarFigureComposite::Pointer pfc = mitk::PlanarFigureComposite::New();
if (operation.empty())
{
result = inputTractogram->ExtractFiberSubset(input1);
}
else if (operation=="NOT")
{
pfc->setOperationType(mitk::PFCOMPOSITION_NOT_OPERATION);
pfc->addPlanarFigure(input1);
result = inputTractogram->ExtractFiberSubset(pfc);
}
else if (operation=="AND" && pf2.IsNotNull())
{
pfc->setOperationType(mitk::PFCOMPOSITION_AND_OPERATION);
pfc->addPlanarFigure(input1);
pfc->addPlanarFigure(input2);
result = inputTractogram->ExtractFiberSubset(pfc);
}
else if (operation=="OR" && pf2.IsNotNull())
{
pfc->setOperationType(mitk::PFCOMPOSITION_OR_OPERATION);
pfc->addPlanarFigure(input1);
pfc->addPlanarFigure(input2);
result = inputTractogram->ExtractFiberSubset(pfc);
}
else
{
MITK_INFO << "Could not process input:";
MITK_INFO << pf1_path;
MITK_INFO << pf2_path;
MITK_INFO << operation;
}
}
else
{
ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New();
mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(pf1_path)->GetData());
mitk::CastToItkImage<ItkUcharImgType>(mitkMaskImage, itkMaskImage);
if (operation=="NOT")
result = inputTractogram->ExtractFiberSubset(itkMaskImage, true, true);
else
result = inputTractogram->ExtractFiberSubset(itkMaskImage, true, false);
}
if (result.IsNotNull())
mitk::IOUtil::SaveBaseData(result, outFib);
else
MITK_INFO << "No valid fiber bundle extracted.";
-
- MITK_INFO << "DONE";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(FiberExtraction);
diff --git a/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp b/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp
index 4cc65437d9..52d81f0705 100644
--- a/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp
+++ b/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp
@@ -1,240 +1,240 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include "MiniAppManager.h"
#include <vector>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <string>
#include <itkImageFileWriter.h>
#include <itkMetaDataObject.h>
#include <itkVectorImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkFiberBundleX.h>
#include "ctkCommandLineParser.h"
#include <boost/lexical_cast.hpp>
#include <mitkCoreObjectFactory.h>
mitk::FiberBundleX::Pointer LoadFib(std::string filename)
{
const std::string s1="", s2="";
std::vector<mitk::BaseData::Pointer> fibInfile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false );
if( fibInfile.empty() )
MITK_INFO << "File " << filename << " could not be read!";
mitk::BaseData::Pointer baseData = fibInfile.at(0);
return dynamic_cast<mitk::FiberBundleX*>(baseData.GetPointer());
}
int FiberProcessing(int argc, char* argv[])
{
+ MITK_INFO << "FiberProcessing";
ctkCommandLineParser parser;
parser.setTitle("Fiber Processing");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib)", us::Any(), false);
parser.addArgument("outFile", "o", ctkCommandLineParser::OutputFile, "Output:", "output fiber bundle (.fib)", us::Any(), false);
parser.addArgument("resample", "r", ctkCommandLineParser::Float, "Resample:", "Resample fiber with the given point distance (in mm)");
parser.addArgument("smooth", "s", ctkCommandLineParser::Float, "Smooth:", "Smooth fiber with the given point distance (in mm)");
parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "Minimum length:", "Minimum fiber length (in mm)");
parser.addArgument("maxLength", "m", ctkCommandLineParser::Float, "Maximum length:", "Maximum fiber length (in mm)");
parser.addArgument("minCurv", "a", ctkCommandLineParser::Float, "Minimum curvature radius:", "Minimum curvature radius (in mm)");
parser.addArgument("mirror", "p", ctkCommandLineParser::Int, "Invert coordinates:", "Invert fiber coordinates XYZ (e.g. 010 to invert y-coordinate of each fiber point)");
parser.addArgument("copyAndJoin", "c", ctkCommandLineParser::Bool, "Copy & Join:", "Create a copy of the input fiber bundle (applied after resample/smooth/minLength/maxLength/minCurv/mirror) and join copy with original (applied after rotate/scale/translate)");
//parser.addArgument("join", "j", ctkCommandLineParser::Bool, "Join the original and copied fiber bundle (applied after rotate/scale/translate)");
parser.addArgument("rotate-x", "rx", ctkCommandLineParser::Float, "Rotate x-axis:", "Rotate around x-axis (if copy is given the copy is rotated, in deg)");
parser.addArgument("rotate-y", "ry", ctkCommandLineParser::Float, "Rotate y-axis:", "Rotate around y-axis (if copy is given the copy is rotated, in deg)");
parser.addArgument("rotate-z", "rz", ctkCommandLineParser::Float, "Rotate z-axis:", "Rotate around z-axis (if copy is given the copy is rotated, in deg)");
parser.addArgument("scale-x", "sx", ctkCommandLineParser::Float, "Scale x-axis:", "Scale in direction of x-axis (if copy is given the copy is scaled)");
parser.addArgument("scale-y", "sy", ctkCommandLineParser::Float, "Scale y-axis:", "Scale in direction of y-axis (if copy is given the copy is scaled)");
parser.addArgument("scale-z", "sz", ctkCommandLineParser::Float, "Scale z-axis", "Scale in direction of z-axis (if copy is given the copy is scaled)");
parser.addArgument("translate-x", "tx", ctkCommandLineParser::Float, "Translate x-axis:", "Translate in direction of x-axis (if copy is given the copy is translated, in mm)");
parser.addArgument("translate-y", "ty", ctkCommandLineParser::Float, "Translate y-axis:", "Translate in direction of y-axis (if copy is given the copy is translated, in mm)");
parser.addArgument("translate-z", "tz", ctkCommandLineParser::Float, "Translate z-axis:", "Translate in direction of z-axis (if copy is given the copy is translated, in mm)");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
float pointDist = -1;
if (parsedArgs.count("resample"))
pointDist = us::any_cast<float>(parsedArgs["resample"]);
float smoothDist = -1;
if (parsedArgs.count("smooth"))
smoothDist = us::any_cast<float>(parsedArgs["smooth"]);
float minFiberLength = -1;
if (parsedArgs.count("minLength"))
minFiberLength = us::any_cast<float>(parsedArgs["minLength"]);
float maxFiberLength = -1;
if (parsedArgs.count("maxLength"))
maxFiberLength = us::any_cast<float>(parsedArgs["maxLength"]);
float curvThres = -1;
if (parsedArgs.count("minCurv"))
curvThres = us::any_cast<float>(parsedArgs["minCurv"]);
int axis = 0;
if (parsedArgs.count("mirror"))
axis = us::any_cast<int>(parsedArgs["mirror"]);
bool copyAndJoin = false;
if(parsedArgs.count("copyAndJoin"))
copyAndJoin = us::any_cast<bool>(parsedArgs["copyAndJoin"]);
float rotateX = 0;
if (parsedArgs.count("rotate-x"))
rotateX = us::any_cast<float>(parsedArgs["rotate-x"]);
float rotateY = 0;
if (parsedArgs.count("rotate-y"))
rotateY = us::any_cast<float>(parsedArgs["rotate-y"]);
float rotateZ = 0;
if (parsedArgs.count("rotate-z"))
rotateZ = us::any_cast<float>(parsedArgs["rotate-z"]);
float scaleX = 0;
if (parsedArgs.count("scale-x"))
scaleX = us::any_cast<float>(parsedArgs["scale-x"]);
float scaleY = 0;
if (parsedArgs.count("scale-y"))
scaleY = us::any_cast<float>(parsedArgs["scale-y"]);
float scaleZ = 0;
if (parsedArgs.count("scale-z"))
scaleZ = us::any_cast<float>(parsedArgs["scale-z"]);
float translateX = 0;
if (parsedArgs.count("translate-x"))
translateX = us::any_cast<float>(parsedArgs["translate-x"]);
float translateY = 0;
if (parsedArgs.count("translate-y"))
translateY = us::any_cast<float>(parsedArgs["translate-y"]);
float translateZ = 0;
if (parsedArgs.count("translate-z"))
translateZ = us::any_cast<float>(parsedArgs["translate-z"]);
string inFileName = us::any_cast<string>(parsedArgs["input"]);
string outFileName = us::any_cast<string>(parsedArgs["outFile"]);
try
{
mitk::FiberBundleX::Pointer fib = LoadFib(inFileName);
if (minFiberLength>0)
fib->RemoveShortFibers(minFiberLength);
if (maxFiberLength>0)
fib->RemoveLongFibers(maxFiberLength);
if (curvThres>0)
fib->ApplyCurvatureThreshold(curvThres, false);
if (pointDist>0)
fib->ResampleFibers(pointDist);
if (smoothDist>0)
fib->DoFiberSmoothing(smoothDist);
if (axis/100==1)
fib->MirrorFibers(0);
if ((axis%100)/10==1)
fib->MirrorFibers(1);
if (axis%10==1)
fib->MirrorFibers(2);
if (copyAndJoin == true)
{
MITK_INFO << "Create copy";
mitk::FiberBundleX::Pointer fibCopy = fib->GetDeepCopy();
if (rotateX > 0 || rotateY > 0 || rotateZ > 0){
MITK_INFO << "Rotate " << rotateX << " " << rotateY << " " << rotateZ;
fibCopy->RotateAroundAxis(rotateX, rotateY, rotateZ);
}
if (translateX > 0 || translateY > 0 || translateZ > 0)
fibCopy->TranslateFibers(translateX, translateY, translateZ);
if (scaleX > 0 || scaleY > 0 || scaleZ > 0)
fibCopy->ScaleFibers(scaleX, scaleY, scaleZ);
MITK_INFO << "Join copy with original";
fib = fib->AddBundle(fibCopy.GetPointer());
} else {
if (rotateX > 0 || rotateY > 0 || rotateZ > 0){
MITK_INFO << "Rotate " << rotateX << " " << rotateY << " " << rotateZ;
fib->RotateAroundAxis(rotateX, rotateY, rotateZ);
}
if (translateX > 0 || translateY > 0 || translateZ > 0){
fib->TranslateFibers(translateX, translateY, translateZ);
}
if (scaleX > 0 || scaleY > 0 || scaleZ > 0)
fib->ScaleFibers(scaleX, scaleY, scaleZ);
}
mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters();
for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it)
{
if ( (*it)->CanWriteBaseDataType(fib.GetPointer()) ) {
- MITK_INFO << "writing " << outFileName;
(*it)->SetFileName( outFileName.c_str() );
(*it)->DoWrite( fib.GetPointer() );
}
}
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(FiberProcessing);
diff --git a/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp b/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp
index 8dfba05cc5..b433de1774 100755
--- a/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp
+++ b/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp
@@ -1,84 +1,85 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <mitkDiffusionImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkIOUtil.h>
#include <mitkNrrdDiffusionImageWriter.h>
#include <mitkFiberBundleX.h>
#include <mitkFiberfoxParameters.h>
#include "ctkCommandLineParser.h"
#include <itkAddArtifactsToDwiImageFilter.h>
#include <itkTractsToDWIImageFilter.h>
#include "boost/property_tree/ptree.hpp"
#include "boost/property_tree/xml_parser.hpp"
#include "boost/foreach.hpp"
/** TODO: Proritype signal komplett speichern oder bild mit speichern. */
/** TODO: Tarball aus images und parametern? */
/** TODO: Artefakte auf bild in miniapp */
namespace mitk
{
int Fiberfox(int argc, char* argv[])
{
+ MITK_INFO << "Fiberfox";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("out", "o", ctkCommandLineParser::OutputFile, "Output root:", "output root", us::Any(), false);
parser.addArgument("parameters", "p", ctkCommandLineParser::InputFile, "Parameter file:", "fiberfox parameter file", us::Any(), false);
parser.addArgument("fiberbundle", "f", ctkCommandLineParser::String, "Fiberbundle:", "", us::Any(), false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string outName = us::any_cast<string>(parsedArgs["out"]);
string paramName = us::any_cast<string>(parsedArgs["parameters"]);
string fibFile = "";
if (parsedArgs.count("fiberbundle"))
fibFile = us::any_cast<string>(parsedArgs["fiberbundle"]);
{
FiberfoxParameters<double> parameters;
parameters.LoadParameters(paramName);
mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(mitk::IOUtil::LoadDataNode(fibFile)->GetData());
itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New();
tractsToDwiFilter->SetParameters(parameters);
tractsToDwiFilter->SetFiberBundle(inputTractogram);
tractsToDwiFilter->Update();
DiffusionImage<short>::Pointer image = DiffusionImage<short>::New();
image->SetVectorImage( tractsToDwiFilter->GetOutput() );
image->SetReferenceBValue( parameters.m_SignalGen.m_Bvalue );
image->SetDirections( parameters.m_SignalGen.GetGradientDirections() );
image->InitializeFromVectorImage();
NrrdDiffusionImageWriter<short>::Pointer writer = NrrdDiffusionImageWriter<short>::New();
writer->SetFileName(outName);
writer->SetInput(image);
writer->Update();
}
return EXIT_SUCCESS;
}
}
RegisterDiffusionMiniApp(Fiberfox);
diff --git a/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp b/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp
index 675b71e0f2..896b1456de 100755
--- a/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp
+++ b/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp
@@ -1,100 +1,96 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <mitkDiffusionImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkIOUtil.h>
#include <mitkNrrdDiffusionImageWriter.h>
#include <mitkFiberBundleX.h>
#include <mitkFiberBundleXWriter.h>
#include "ctkCommandLineParser.h"
#include "ctkCommandLineParser.cpp"
using namespace mitk;
int FileFormatConverter(int argc, char* argv[])
{
ctkCommandLineParser parser;
parser.setTitle("Format Converter");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("in", "i", ctkCommandLineParser::InputFile, "Input:", "input file", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::OutputFile, "Output:", "output file", us::Any(), false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
// mandatory arguments
string inName = us::any_cast<string>(parsedArgs["in"]);
string outName = us::any_cast<string>(parsedArgs["out"]);
try
{
- MITK_INFO << "Loading " << inName;
const std::string s1="", s2="";
std::vector<BaseData::Pointer> infile = BaseDataIO::LoadBaseDataFromFile( inName, s1, s2, false );
mitk::BaseData::Pointer baseData = infile.at(0);
if ( dynamic_cast<DiffusionImage<short>*>(baseData.GetPointer()) )
{
- MITK_INFO << "Writing " << outName;
DiffusionImage<short>::Pointer dwi = dynamic_cast<DiffusionImage<short>*>(baseData.GetPointer());
NrrdDiffusionImageWriter<short>::Pointer writer = NrrdDiffusionImageWriter<short>::New();
writer->SetFileName(outName);
writer->SetInput(dwi);
writer->Update();
}
else if ( dynamic_cast<Image*>(baseData.GetPointer()) )
{
Image::Pointer image = dynamic_cast<Image*>(baseData.GetPointer());
mitk::IOUtil::SaveImage(image, outName);
}
else if ( dynamic_cast<FiberBundleX*>(baseData.GetPointer()) )
{
- MITK_INFO << "Writing " << outName;
FiberBundleXWriter::Pointer fibWriter = FiberBundleXWriter::New();
fibWriter->SetFileName(outName.c_str());
fibWriter->DoWrite( dynamic_cast<FiberBundleX*>(baseData.GetPointer()) );
}
else
MITK_INFO << "File type currently not supported!";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
- MITK_INFO << "DONE";
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(FileFormatConverter);
diff --git a/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp b/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp
index 1c1de05216..65cc3c6046 100755
--- a/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp
+++ b/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp
@@ -1,251 +1,252 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <mitkQBallImage.h>
#include <mitkTensorImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkFiberBundleX.h>
#include <itkGibbsTrackingFilter.h>
#include <itkDiffusionTensor3D.h>
#include <itkShCoefficientImageImporter.h>
#include <mitkImageToItk.h>
#include <mitkIOUtil.h>
#include "ctkCommandLineParser.h"
#include <boost/algorithm/string.hpp>
#include <itkFlipImageFilter.h>
#include <mitkCoreObjectFactory.h>
template<int shOrder>
typename itk::ShCoefficientImageImporter< float, shOrder >::QballImageType::Pointer TemplatedConvertShCoeffs(mitk::Image* mitkImg, int toolkit, bool noFlip = false)
{
typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType;
typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType;
CasterType::Pointer caster = CasterType::New();
caster->SetInput(mitkImg);
caster->Update();
itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput();
typename FilterType::Pointer filter = FilterType::New();
if (noFlip)
{
filter->SetInputImage(itkImage);
}
else
{
MITK_INFO << "Flipping image";
itk::FixedArray<bool, 4> flipAxes;
flipAxes[0] = true;
flipAxes[1] = true;
flipAxes[2] = false;
flipAxes[3] = false;
itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New();
flipper->SetInput(itkImage);
flipper->SetFlipAxes(flipAxes);
flipper->Update();
itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput();
itk::Matrix< double,4,4 > m = itkImage->GetDirection(); m[0][0] *= -1; m[1][1] *= -1;
flipped->SetDirection(m);
itk::Point< float, 4 > o = itkImage->GetOrigin();
o[0] -= (flipped->GetLargestPossibleRegion().GetSize(0)-1);
o[1] -= (flipped->GetLargestPossibleRegion().GetSize(1)-1);
flipped->SetOrigin(o);
filter->SetInputImage(flipped);
}
switch (toolkit)
{
case 0:
filter->SetToolkit(FilterType::FSL);
break;
case 1:
filter->SetToolkit(FilterType::MRTRIX);
break;
default:
filter->SetToolkit(FilterType::FSL);
}
filter->GenerateData();
return filter->GetQballImage();
}
int GibbsTracking(int argc, char* argv[])
{
+ MITK_INFO << "GibbsTracking";
ctkCommandLineParser parser;
parser.setTitle("Gibbs Tracking");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input:", "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", us::Any(), false);
parser.addArgument("parameters", "p", ctkCommandLineParser::InputFile, "Parameters:", "parameter file (.gtp)", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::InputFile, "Mask:", "binary mask image");
parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "SH coefficient:", "sh coefficient convention (FSL, MRtrix)", string("FSL"), true);
parser.addArgument("outFile", "o", ctkCommandLineParser::OutputFile, "Output:", "output fiber bundle (.fib)", us::Any(), false);
parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "No flip:", "do not flip input image to match MITK coordinate convention");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string inFileName = us::any_cast<string>(parsedArgs["input"]);
string paramFileName = us::any_cast<string>(parsedArgs["parameters"]);
string outFileName = us::any_cast<string>(parsedArgs["outFile"]);
bool noFlip = false;
if (parsedArgs.count("noFlip"))
noFlip = us::any_cast<bool>(parsedArgs["noFlip"]);
try
{
// instantiate gibbs tracker
typedef itk::Vector<float, QBALL_ODFSIZE> OdfVectorType;
typedef itk::Image<OdfVectorType,3> ItkQballImageType;
typedef itk::GibbsTrackingFilter<ItkQballImageType> GibbsTrackingFilterType;
GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New();
// load input image
const std::string s1="", s2="";
std::vector<mitk::BaseData::Pointer> infile = mitk::BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false );
mitk::Image::Pointer mitkImage = dynamic_cast<mitk::Image*>(infile.at(0).GetPointer());
// try to cast to qball image
if( boost::algorithm::ends_with(inFileName, ".qbi") )
{
MITK_INFO << "Loading qball image ...";
mitk::QBallImage::Pointer mitkQballImage = dynamic_cast<mitk::QBallImage*>(infile.at(0).GetPointer());
ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New();
mitk::CastToItkImage(mitkQballImage, itk_qbi);
gibbsTracker->SetQBallImage(itk_qbi.GetPointer());
}
else if( boost::algorithm::ends_with(inFileName, ".dti") )
{
MITK_INFO << "Loading tensor image ...";
typedef itk::Image< itk::DiffusionTensor3D<float>, 3 > ItkTensorImage;
mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast<mitk::TensorImage*>(infile.at(0).GetPointer());
ItkTensorImage::Pointer itk_dti = ItkTensorImage::New();
mitk::CastToItkImage(mitkTensorImage, itk_dti);
gibbsTracker->SetTensorImage(itk_dti);
}
else if ( boost::algorithm::ends_with(inFileName, ".nii") )
{
MITK_INFO << "Loading sh-coefficient image ...";
int nrCoeffs = mitkImage->GetLargestPossibleRegion().GetSize()[3];
int c=3, d=2-2*nrCoeffs;
double D = c*c-4*d;
int shOrder;
if (D>0)
{
shOrder = (-c+sqrt(D))/2.0;
if (shOrder<0)
shOrder = (-c-sqrt(D))/2.0;
}
else if (D==0)
shOrder = -c/2.0;
MITK_INFO << "using SH-order " << shOrder;
int toolkitConvention = 0;
if (parsedArgs.count("shConvention"))
{
string convention = us::any_cast<string>(parsedArgs["shConvention"]).c_str();
if ( boost::algorithm::equals(convention, "MRtrix") )
{
toolkitConvention = 1;
MITK_INFO << "Using MRtrix style sh-coefficient convention";
}
else
MITK_INFO << "Using FSL style sh-coefficient convention";
}
else
MITK_INFO << "Using FSL style sh-coefficient convention";
switch (shOrder)
{
case 4:
gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<4>(mitkImage, toolkitConvention, noFlip));
break;
case 6:
gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<6>(mitkImage, toolkitConvention, noFlip));
break;
case 8:
gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<8>(mitkImage, toolkitConvention, noFlip));
break;
case 10:
gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<10>(mitkImage, toolkitConvention, noFlip));
break;
case 12:
gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<12>(mitkImage, toolkitConvention, noFlip));
break;
default:
MITK_INFO << "SH-order " << shOrder << " not supported";
}
}
else
return EXIT_FAILURE;
// global tracking
if (parsedArgs.count("mask"))
{
typedef itk::Image<float,3> MaskImgType;
mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(us::any_cast<string>(parsedArgs["mask"]));
MaskImgType::Pointer itk_mask = MaskImgType::New();
mitk::CastToItkImage(mitkMaskImage, itk_mask);
gibbsTracker->SetMaskImage(itk_mask);
}
gibbsTracker->SetDuplicateImage(false);
gibbsTracker->SetLoadParameterFile( paramFileName );
// gibbsTracker->SetLutPath( "" );
gibbsTracker->Update();
mitk::FiberBundleX::Pointer mitkFiberBundle = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle());
mitkFiberBundle->SetReferenceImage(mitkImage);
mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters();
for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it)
{
if ( (*it)->CanWriteBaseDataType(mitkFiberBundle.GetPointer()) ) {
(*it)->SetFileName( outFileName.c_str() );
(*it)->DoWrite( mitkFiberBundle.GetPointer() );
}
}
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(GibbsTracking);
diff --git a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp
index 9c5a6e5871..37d9aa0800 100755
--- a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp
+++ b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp
@@ -1,314 +1,309 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <itkEvaluateDirectionImagesFilter.h>
#include <metaCommand.h>
#include "ctkCommandLineParser.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 LocalDirectionalFiberPlausibility(int argc, char* argv[])
{
+ MITK_INFO << "LocalDirectionalFiberPlausibility";
ctkCommandLineParser 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", ctkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "Reference images:", "reference direction images", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::StringList, "Masks:", "mask images");
parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true);
parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results");
parser.addArgument("ignore", "n", ctkCommandLineParser::Bool, "Ignore:", "don't increase error for missing or too many directions");
parser.addArgument("fileID", "id", ctkCommandLineParser::String, "ID:", "optional ID field");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
ctkCommandLineParser::StringContainerType referenceImages = us::any_cast<ctkCommandLineParser::StringContainerType>(parsedArgs["reference"]);
ctkCommandLineParser::StringContainerType maskImages;
if (parsedArgs.count("mask"))
maskImages = us::any_cast<ctkCommandLineParser::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"]);
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;
if (parsedArgs.count("ignore"))
ignore = us::any_cast<bool>(parsedArgs["ignore"]);
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::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(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(...){ MITK_INFO << "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->Update();
ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer();
if (verbose)
{
// write vector field
mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle();
mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters();
for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it)
{
if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) {
string outfilename = outRoot;
outfilename.append("_VECTOR_FIELD.fib");
(*it)->SetFileName( outfilename.c_str() );
(*it)->DoWrite( directions.GetPointer() );
}
}
// 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");
- MITK_INFO << "writing " << outfilename;
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");
- MITK_INFO << "writing " << outfilename;
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->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");
- MITK_INFO << "writing " << outfilename;
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::GetFilenameWithoutExtension(itksys::SystemTools::GetFilenameName(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->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");
- MITK_INFO << "writing " << outfilename;
writer->SetFileName(outfilename.c_str());
writer->SetInput(angularErrorImage);
writer->Update();
}
string sens = itksys::SystemTools::GetFilenameWithoutExtension(itksys::SystemTools::GetFilenameName(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();
-
- MITK_INFO << "DONE";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(LocalDirectionalFiberPlausibility);
diff --git a/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp b/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp
index 0368b6b6a0..3a81eabccf 100644
--- a/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp
+++ b/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp
@@ -1,230 +1,229 @@
/*===================================================================
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 "MiniAppManager.h"
#include <vector>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <string>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkExceptionObject.h>
#include <itkImageFileWriter.h>
#include <itkMetaDataObject.h>
#include <itkVectorImage.h>
#include <itkResampleImageFilter.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkDiffusionImage.h>
#include <mitkQBallImage.h>
#include <mitkBaseData.h>
#include <mitkFiberBundleX.h>
#include "ctkCommandLineParser.h"
#include <boost/lexical_cast.hpp>
#include <itkRadialMultishellToSingleshellImageFilter.h>
#include <itkADCAverageFunctor.h>
#include <itkBiExpFitFunctor.h>
#include <itkKurtosisFitFunctor.h>
#include <mitkNrrdDiffusionImageWriter.h>
#include <itkDwiGradientLengthCorrectionFilter.h>
int MultishellMethods(int argc, char* argv[])
{
+ MITK_INFO << "MultishellMethods";
ctkCommandLineParser parser;
parser.setTitle("Multishell Methods");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("in", "i", ctkCommandLineParser::InputFile, "Input:", "input file", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::OutputFile, "Output:", "output file", us::Any(), false);
parser.addArgument("adc", "D", ctkCommandLineParser::Bool, "ADC:", "ADC Average", us::Any(), false);
parser.addArgument("akc", "K", ctkCommandLineParser::Bool, "Kurtosis fit:", "Kurtosis Fit", us::Any(), false);
parser.addArgument("biexp", "B", ctkCommandLineParser::Bool, "BiExp fit:", "BiExp fit", us::Any(), false);
parser.addArgument("targetbvalue", "b", ctkCommandLineParser::String, "b Value:", "target bValue (mean, min, max)", us::Any(), false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
// mandatory arguments
string inName = us::any_cast<string>(parsedArgs["in"]);
string outName = us::any_cast<string>(parsedArgs["out"]);
bool applyADC = us::any_cast<bool>(parsedArgs["adc"]);
bool applyAKC = us::any_cast<bool>(parsedArgs["akc"]);
bool applyBiExp = us::any_cast<bool>(parsedArgs["biexp"]);
string targetType = us::any_cast<string>(parsedArgs["targetbvalue"]);
try
{
MITK_INFO << "Loading " << inName;
const std::string s1="", s2="";
std::vector<mitk::BaseData::Pointer> infile = mitk::BaseDataIO::LoadBaseDataFromFile( inName, s1, s2, false );
mitk::BaseData::Pointer baseData = infile.at(0);
if ( dynamic_cast<mitk::DiffusionImage<short>*>(baseData.GetPointer()) )
{
- MITK_INFO << "Writing " << outName;
mitk::DiffusionImage<short>::Pointer dwi = dynamic_cast<mitk::DiffusionImage<short>*>(baseData.GetPointer());
typedef itk::RadialMultishellToSingleshellImageFilter<short, short> FilterType;
typedef itk::DwiGradientLengthCorrectionFilter CorrectionFilterType;
CorrectionFilterType::Pointer roundfilter = CorrectionFilterType::New();
roundfilter->SetRoundingValue( 1000 );
roundfilter->SetReferenceBValue(dwi->GetReferenceBValue());
roundfilter->SetReferenceGradientDirectionContainer(dwi->GetDirections());
roundfilter->Update();
dwi->SetReferenceBValue( roundfilter->GetNewBValue() );
dwi->SetDirections( roundfilter->GetOutputGradientDirectionContainer());
// filter input parameter
const mitk::DiffusionImage<short>::BValueMap
&originalShellMap = dwi->GetBValueMap();
const mitk::DiffusionImage<short>::ImageType
*vectorImage = dwi->GetVectorImage();
const mitk::DiffusionImage<short>::GradientDirectionContainerType::Pointer
gradientContainer = dwi->GetDirections();
const unsigned int
&bValue = dwi->GetReferenceBValue();
// filter call
vnl_vector<double> bValueList(originalShellMap.size()-1);
double targetBValue = bValueList.mean();
mitk::DiffusionImage<short>::BValueMap::const_iterator it = originalShellMap.begin();
++it; int i = 0 ;
for(; it != originalShellMap.end(); ++it)
bValueList.put(i++,it->first);
if( targetType == "mean" )
targetBValue = bValueList.mean();
else if( targetType == "min" )
targetBValue = bValueList.min_value();
else if( targetType == "max" )
targetBValue = bValueList.max_value();
if(applyADC)
{
FilterType::Pointer filter = FilterType::New();
filter->SetInput(vectorImage);
filter->SetOriginalGradientDirections(gradientContainer);
filter->SetOriginalBValueMap(originalShellMap);
filter->SetOriginalBValue(bValue);
itk::ADCAverageFunctor::Pointer functor = itk::ADCAverageFunctor::New();
functor->setListOfBValues(bValueList);
functor->setTargetBValue(targetBValue);
filter->SetFunctor(functor);
filter->Update();
// create new DWI image
mitk::DiffusionImage<short>::Pointer outImage = mitk::DiffusionImage<short>::New();
outImage->SetVectorImage( filter->GetOutput() );
outImage->SetReferenceBValue( targetBValue );
outImage->SetDirections( filter->GetTargetGradientDirections() );
outImage->InitializeFromVectorImage();
mitk::NrrdDiffusionImageWriter<short>::Pointer writer = mitk::NrrdDiffusionImageWriter<short>::New();
writer->SetFileName((string(outName) + "_ADC.dwi"));
writer->SetInput(outImage);
writer->Update();
}
if(applyAKC)
{
FilterType::Pointer filter = FilterType::New();
filter->SetInput(vectorImage);
filter->SetOriginalGradientDirections(gradientContainer);
filter->SetOriginalBValueMap(originalShellMap);
filter->SetOriginalBValue(bValue);
itk::KurtosisFitFunctor::Pointer functor = itk::KurtosisFitFunctor::New();
functor->setListOfBValues(bValueList);
functor->setTargetBValue(targetBValue);
filter->SetFunctor(functor);
filter->Update();
// create new DWI image
mitk::DiffusionImage<short>::Pointer outImage = mitk::DiffusionImage<short>::New();
outImage->SetVectorImage( filter->GetOutput() );
outImage->SetReferenceBValue( targetBValue );
outImage->SetDirections( filter->GetTargetGradientDirections() );
outImage->InitializeFromVectorImage();
mitk::NrrdDiffusionImageWriter<short>::Pointer writer = mitk::NrrdDiffusionImageWriter<short>::New();
writer->SetFileName((string(outName) + "_AKC.dwi"));
writer->SetInput(outImage);
writer->Update();
}
if(applyBiExp)
{
FilterType::Pointer filter = FilterType::New();
filter->SetInput(vectorImage);
filter->SetOriginalGradientDirections(gradientContainer);
filter->SetOriginalBValueMap(originalShellMap);
filter->SetOriginalBValue(bValue);
itk::BiExpFitFunctor::Pointer functor = itk::BiExpFitFunctor::New();
functor->setListOfBValues(bValueList);
functor->setTargetBValue(targetBValue);
filter->SetFunctor(functor);
filter->Update();
// create new DWI image
mitk::DiffusionImage<short>::Pointer outImage = mitk::DiffusionImage<short>::New();
outImage->SetVectorImage( filter->GetOutput() );
outImage->SetReferenceBValue( targetBValue );
outImage->SetDirections( filter->GetTargetGradientDirections() );
outImage->InitializeFromVectorImage();
mitk::NrrdDiffusionImageWriter<short>::Pointer writer = mitk::NrrdDiffusionImageWriter<short>::New();
writer->SetFileName((string(outName) + "_BiExp.dwi"));
writer->SetInput(outImage);
writer->Update();
}
}
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
- MITK_INFO << "DONE";
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(MultishellMethods);
diff --git a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp
index 99e1acb9af..6104af810f 100755
--- a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp
@@ -1,382 +1,380 @@
/*===================================================================
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 "MiniAppManager.h"
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkFiniteDiffOdfMaximaExtractionFilter.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkDiffusionImage.h>
#include <mitkQBallImage.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <mitkTensorImage.h>
#include <mitkCoreObjectFactory.h>
#include "ctkCommandLineParser.h"
#include <mitkFiberBundleXWriter.h>
#include <itkShCoefficientImageImporter.h>
#include <itkFlipImageFilter.h>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
mitk::Image::Pointer LoadData(std::string filename)
{
if( filename.empty() )
return NULL;
const std::string s1="", s2="";
std::vector<mitk::BaseData::Pointer> infile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false );
if( infile.empty() )
{
MITK_INFO << "File " << filename << " could not be read!";
return NULL;
}
mitk::BaseData::Pointer baseData = infile.at(0);
return dynamic_cast<mitk::Image*>(baseData.GetPointer());
}
template<int shOrder>
int StartPeakExtraction(int argc, char* argv[])
{
+ MITK_INFO << "StartPeakExtraction";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("image", "i", ctkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false);
parser.addArgument("outroot", "o", ctkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::InputFile, "Mask", "mask image");
parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true);
parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true);
parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true);
parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true);
parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true);
parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention");
parser.setCategory("Preprocessing Tools");
parser.setTitle("Peak Extraction");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
// mandatory arguments
string imageName = us::any_cast<string>(parsedArgs["image"]);
string outRoot = us::any_cast<string>(parsedArgs["outroot"]);
// optional arguments
string maskImageName("");
if (parsedArgs.count("mask"))
maskImageName = us::any_cast<string>(parsedArgs["mask"]);
int normalization = 1;
if (parsedArgs.count("normalization"))
normalization = us::any_cast<int>(parsedArgs["normalization"]);
int numPeaks = 2;
if (parsedArgs.count("numpeaks"))
numPeaks = us::any_cast<int>(parsedArgs["numpeaks"]);
float peakThres = 0.4;
if (parsedArgs.count("peakthres"))
peakThres = us::any_cast<float>(parsedArgs["peakthres"]);
float absPeakThres = 0.06;
if (parsedArgs.count("abspeakthres"))
absPeakThres = us::any_cast<float>(parsedArgs["abspeakthres"]);
bool noFlip = false;
if (parsedArgs.count("noFlip"))
noFlip = us::any_cast<bool>(parsedArgs["noFlip"]);
MITK_INFO << "image: " << imageName;
MITK_INFO << "outroot: " << outRoot;
if (!maskImageName.empty())
MITK_INFO << "mask: " << maskImageName;
else
MITK_INFO << "no mask image selected";
MITK_INFO << "numpeaks: " << numPeaks;
MITK_INFO << "peakthres: " << peakThres;
MITK_INFO << "abspeakthres: " << absPeakThres;
MITK_INFO << "shOrder: " << shOrder;
try
{
mitk::Image::Pointer image = LoadData(imageName);
mitk::Image::Pointer mask = LoadData(maskImageName);
typedef itk::Image<unsigned char, 3> ItkUcharImgType;
typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType;
typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New();
int toolkitConvention = 0;
if (parsedArgs.count("shConvention"))
{
string convention = us::any_cast<string>(parsedArgs["shConvention"]).c_str();
if ( boost::algorithm::equals(convention, "FSL") )
{
toolkitConvention = 1;
MITK_INFO << "Using FSL SH-basis";
}
else if ( boost::algorithm::equals(convention, "MRtrix") )
{
toolkitConvention = 2;
MITK_INFO << "Using MRtrix SH-basis";
}
else
MITK_INFO << "Using MITK SH-basis";
}
else
MITK_INFO << "Using MITK SH-basis";
ItkUcharImgType::Pointer itkMaskImage = NULL;
if (mask.IsNotNull())
{
try{
itkMaskImage = ItkUcharImgType::New();
mitk::CastToItkImage(mask, itkMaskImage);
filter->SetMaskImage(itkMaskImage);
}
catch(...)
{
}
}
if (toolkitConvention>0)
{
MITK_INFO << "Converting coefficient image to MITK format";
typedef itk::ShCoefficientImageImporter< float, shOrder > ConverterType;
typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType;
CasterType::Pointer caster = CasterType::New();
caster->SetInput(image);
caster->Update();
itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput();
typename ConverterType::Pointer converter = ConverterType::New();
if (noFlip)
{
converter->SetInputImage(itkImage);
}
else
{
MITK_INFO << "Flipping image";
itk::FixedArray<bool, 4> flipAxes;
flipAxes[0] = true;
flipAxes[1] = true;
flipAxes[2] = false;
flipAxes[3] = false;
itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New();
flipper->SetInput(itkImage);
flipper->SetFlipAxes(flipAxes);
flipper->Update();
itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput();
itk::Matrix< double,4,4 > m = itkImage->GetDirection(); m[0][0] *= -1; m[1][1] *= -1;
flipped->SetDirection(m);
itk::Point< float, 4 > o = itkImage->GetOrigin();
o[0] -= (flipped->GetLargestPossibleRegion().GetSize(0)-1);
o[1] -= (flipped->GetLargestPossibleRegion().GetSize(1)-1);
flipped->SetOrigin(o);
converter->SetInputImage(flipped);
}
MITK_INFO << "Starting conversion";
switch (toolkitConvention)
{
case 1:
converter->SetToolkit(ConverterType::FSL);
filter->SetToolkit(MaximaExtractionFilterType::FSL);
break;
case 2:
converter->SetToolkit(ConverterType::MRTRIX);
filter->SetToolkit(MaximaExtractionFilterType::MRTRIX);
break;
default:
converter->SetToolkit(ConverterType::FSL);
filter->SetToolkit(MaximaExtractionFilterType::FSL);
break;
}
converter->GenerateData();
filter->SetInput(converter->GetCoefficientImage());
}
else
{
try{
typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType;
typename CasterType::Pointer caster = CasterType::New();
caster->SetInput(image);
caster->Update();
filter->SetInput(caster->GetOutput());
}
catch(...)
{
MITK_INFO << "wrong image type";
return EXIT_FAILURE;
}
}
filter->SetMaxNumPeaks(numPeaks);
filter->SetPeakThreshold(peakThres);
filter->SetAbsolutePeakThreshold(absPeakThres);
filter->SetAngularThreshold(1);
switch (normalization)
{
case 0:
filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM);
break;
case 1:
filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM);
break;
case 2:
filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM);
break;
}
MITK_INFO << "Starting extraction";
filter->Update();
// write direction images
{
typedef typename MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer;
typename ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer();
for (unsigned int i=0; i<container->Size(); i++)
{
typename MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i);
if (itkMaskImage.IsNotNull())
{
itkImg->SetDirection(itkMaskImage->GetDirection());
itkImg->SetOrigin(itkMaskImage->GetOrigin());
}
string outfilename = outRoot;
outfilename.append("_DIRECTION_");
outfilename.append(boost::lexical_cast<string>(i));
outfilename.append(".nrrd");
- MITK_INFO << "writing " << outfilename;
typedef itk::ImageFileWriter< typename MaximaExtractionFilterType::ItkDirectionImage > WriterType;
typename WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outfilename);
writer->SetInput(itkImg);
writer->Update();
}
}
// write num directions image
{
ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage();
if (itkMaskImage.IsNotNull())
{
numDirImage->SetDirection(itkMaskImage->GetDirection());
numDirImage->SetOrigin(itkMaskImage->GetOrigin());
}
string outfilename = outRoot.c_str();
outfilename.append("_NUM_DIRECTIONS.nrrd");
- MITK_INFO << "writing " << outfilename;
typedef itk::ImageFileWriter< ItkUcharImgType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outfilename);
writer->SetInput(numDirImage);
writer->Update();
}
// write vector field
{
mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle();
string outfilename = outRoot.c_str();
outfilename.append("_VECTOR_FIELD.fib");
mitk::FiberBundleXWriter::Pointer fibWriter = mitk::FiberBundleXWriter::New();
fibWriter->SetFileName(outfilename.c_str());
fibWriter->DoWrite(directions.GetPointer());
}
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
- MITK_INFO << "DONE";
return EXIT_SUCCESS;
}
int PeakExtraction(int argc, char* argv[])
{
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("image", "i", ctkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false);
parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order");
parser.addArgument("outroot", "o", ctkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::InputFile, "Mask", "mask image");
parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true);
parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true);
parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true);
parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true);
parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true);
parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention");
parser.setCategory("Preprocessing Tools");
parser.setTitle("Peak Extraction");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
int shOrder = -1;
if (parsedArgs.count("shOrder"))
shOrder = us::any_cast<int>(parsedArgs["shOrder"]);
switch (shOrder)
{
case 4:
return StartPeakExtraction<4>(argc, argv);
case 6:
return StartPeakExtraction<6>(argc, argv);
case 8:
return StartPeakExtraction<8>(argc, argv);
case 10:
return StartPeakExtraction<10>(argc, argv);
case 12:
return StartPeakExtraction<12>(argc, argv);
}
return EXIT_FAILURE;
}
RegisterDiffusionMiniApp(PeakExtraction);
diff --git a/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp b/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp
index b3747772f2..9b5d4bbe10 100755
--- a/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp
+++ b/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp
@@ -1,210 +1,208 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <itkEvaluateDirectionImagesFilter.h>
#include <metaCommand.h>
#include "ctkCommandLineParser.h"
#include <itkTractsToVectorImageFilter.h>
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
#define _USE_MATH_DEFINES
#include <math.h>
int PeaksAngularError(int argc, char* argv[])
{
+ MITK_INFO << "PeaksAngularError";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("test", "t", ctkCommandLineParser::StringList, "Test images", "test direction images", us::Any(), false);
parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "Reference images", "reference direction images", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::InputFile, "Mask", "mask image");
parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "Verbose", "output optional and intermediate calculation results");
parser.addArgument("ignore", "i", ctkCommandLineParser::Bool, "Ignore", "don't increase error for missing or too many directions");
parser.setCategory("Preprocessing Tools");
parser.setTitle("Peaks Angular Error");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
ctkCommandLineParser::StringContainerType testImages = us::any_cast<ctkCommandLineParser::StringContainerType>(parsedArgs["test"]);
ctkCommandLineParser::StringContainerType referenceImages = us::any_cast<ctkCommandLineParser::StringContainerType>(parsedArgs["reference"]);
string maskImage("");
if (parsedArgs.count("mask"))
maskImage = us::any_cast<string>(parsedArgs["mask"]);
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;
if (parsedArgs.count("ignore"))
ignore = us::any_cast<bool>(parsedArgs["ignore"]);
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;
ItkDirectionImageContainerType::Pointer directionImageContainer = ItkDirectionImageContainerType::New();
for (unsigned int i=0; i<testImages.size(); i++)
{
try
{
mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(testImages.at(i))->GetData());
typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType;
CasterType::Pointer caster = CasterType::New();
caster->SetInput(img);
caster->Update();
ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput();
directionImageContainer->InsertElement(directionImageContainer->Size(),itkImg);
}
catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); }
}
// 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(...){ MITK_INFO << "could not load: " << referenceImages.at(i); }
}
// load/create mask image
ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New();
if (maskImage.compare("")==0)
{
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);
}
else
{
mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(maskImage)->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->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");
- MITK_INFO << "writing " << outfilename;
writer->SetFileName(outfilename.c_str());
writer->SetInput(angularErrorImage);
writer->Update();
}
string logFile = outRoot;
logFile.append("_ANGULAR_ERROR.csv");
ofstream file;
file.open (logFile.c_str());
string sens = "Mean:";
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMeanAngularError()));
sens.append(";\n");
sens.append("Median:");
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMedianAngularError()));
sens.append(";\n");
sens.append("Maximum:");
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMaxAngularError()));
sens.append(";\n");
sens.append("Minimum:");
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMinAngularError()));
sens.append(";\n");
sens.append("STDEV:");
sens.append(",");
sens.append(boost::lexical_cast<string>(std::sqrt(evaluationFilter->GetVarAngularError())));
sens.append(";\n");
file << sens;
file.close();
-
- MITK_INFO << "DONE";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(PeaksAngularError);
diff --git a/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp b/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp
index d0531a20ba..68392158d1 100644
--- a/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp
@@ -1,244 +1,243 @@
/*===================================================================
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 "MiniAppManager.h"
#include "mitkBaseDataIOFactory.h"
#include <mitkCoreObjectFactory.h>
#include "mitkDiffusionImage.h"
#include "itkAnalyticalDiffusionQballReconstructionImageFilter.h"
#include <boost/lexical_cast.hpp>
#include "ctkCommandLineParser.h"
#include <mitkIOUtil.h>
#include <itksys/SystemTools.hxx>
using namespace mitk;
/**
* Perform Q-ball reconstruction using a spherical harmonics basis
*/
int QballReconstruction(int argc, char* argv[])
{
+ MITK_INFO << "QballReconstruction";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false);
parser.addArgument("outFile", "o", ctkCommandLineParser::OutputFile, "Output file", "output file", us::Any(), false);
parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order", 4, true);
parser.addArgument("b0Threshold", "t", ctkCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0, true);
parser.addArgument("lambda", "r", ctkCommandLineParser::Float, "Lambda", "ragularization factor lambda", 0.006, true);
parser.addArgument("csa", "csa", ctkCommandLineParser::Bool, "Constant solid angle consideration", "use constant solid angle consideration");
parser.addArgument("outputCoeffs", "shc", ctkCommandLineParser::Bool, "Output coefficients", "output file containing the SH coefficients");
parser.addArgument("mrtrix", "mb", ctkCommandLineParser::Bool, "MRtrix", "use MRtrix compatible spherical harmonics definition");
parser.setCategory("Preprocessing Tools");
parser.setTitle("Qball Reconstruction");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
std::string inFileName = us::any_cast<string>(parsedArgs["input"]);
std::string outfilename = us::any_cast<string>(parsedArgs["outFile"]);
outfilename = itksys::SystemTools::GetFilenamePath(outfilename)+"/"+itksys::SystemTools::GetFilenameWithoutExtension(outfilename);
int threshold = 0;
if (parsedArgs.count("b0Threshold"))
threshold = us::any_cast<int>(parsedArgs["b0Threshold"]);
int shOrder = 4;
if (parsedArgs.count("shOrder"))
shOrder = us::any_cast<int>(parsedArgs["shOrder"]);
float lambda = 0.006;
if (parsedArgs.count("lambda"))
lambda = us::any_cast<float>(parsedArgs["lambda"]);
int normalization = 0;
if (parsedArgs.count("csa") && us::any_cast<bool>(parsedArgs["csa"]))
normalization = 6;
bool outCoeffs = false;
if (parsedArgs.count("outputCoeffs"))
outCoeffs = us::any_cast<bool>(parsedArgs["outputCoeffs"]);
bool mrTrix = false;
if (parsedArgs.count("mrtrix"))
mrTrix = us::any_cast<bool>(parsedArgs["mrtrix"]);
try
{
- MITK_INFO << "Loading image ...";
const std::string s1="", s2="";
std::vector<BaseData::Pointer> infile = BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false );
DiffusionImage<short>::Pointer dwi = dynamic_cast<DiffusionImage<short>*>(infile.at(0).GetPointer());
dwi->AverageRedundantGradients(0.001);
mitk::QBallImage::Pointer image = mitk::QBallImage::New();
mitk::Image::Pointer coeffsImage = mitk::Image::New();
MITK_INFO << "SH order: " << shOrder;
MITK_INFO << "lambda: " << lambda;
MITK_INFO << "B0 threshold: " << threshold;
switch ( shOrder )
{
case 4:
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold( threshold );
filter->SetLambda(lambda);
filter->SetUseMrtrixBasis(mrTrix);
if (normalization==0)
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
else
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
break;
}
case 6:
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,6,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold( threshold );
filter->SetLambda(lambda);
filter->SetUseMrtrixBasis(mrTrix);
if (normalization==0)
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
else
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
break;
}
case 8:
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,8,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold( threshold );
filter->SetLambda(lambda);
filter->SetUseMrtrixBasis(mrTrix);
if (normalization==0)
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
else
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
break;
}
case 10:
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,10,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold( threshold );
filter->SetLambda(lambda);
filter->SetUseMrtrixBasis(mrTrix);
if (normalization==0)
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
else
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
break;
}
case 12:
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,12,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold( threshold );
filter->SetLambda(lambda);
if (normalization==0)
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
else
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
break;
}
default:
{
MITK_INFO << "Supplied SH order not supported. Using default order of 4.";
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold( threshold );
filter->SetLambda(lambda);
filter->SetUseMrtrixBasis(mrTrix);
if (normalization==0)
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
else
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() );
coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() );
}
}
std::string coeffout = outfilename;
coeffout += "_shcoeffs.nrrd";
outfilename += ".qbi";
- MITK_INFO << "writing image " << outfilename;
mitk::IOUtil::SaveBaseData(image, outfilename);
if (outCoeffs)
mitk::IOUtil::SaveImage(coeffsImage, coeffout);
}
catch ( itk::ExceptionObject &err)
{
MITK_INFO << "Exception: " << err;
}
catch ( std::exception err)
{
MITK_INFO << "Exception: " << err.what();
}
catch ( ... )
{
MITK_INFO << "Exception!";
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(QballReconstruction);
diff --git a/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp b/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp
index 6d78c197e7..1e4887d1a7 100755
--- a/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp
+++ b/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp
@@ -1,187 +1,187 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkImageCast.h>
#include <mitkTensorImage.h>
#include <mitkIOUtil.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkFiberBundleX.h>
#include <itkStreamlineTrackingFilter.h>
#include <itkDiffusionTensor3D.h>
#include "ctkCommandLineParser.h"
#include <mitkCoreObjectFactory.h>
#include <mitkFiberBundleXWriter.h>
int StreamlineTracking(int argc, char* argv[])
{
+ MITK_INFO << "StreamlineTracking";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::StringList, "Input image", "input tensor image (.dti)", us::Any(), false);
parser.addArgument("seed", "si", ctkCommandLineParser::InputFile, "Seed image", "binary seed image", us::Any(), true);
parser.addArgument("mask", "mi", ctkCommandLineParser::InputFile, "Mask", "binary mask image", us::Any(), true);
parser.addArgument("faImage", "fai", ctkCommandLineParser::InputFile, "FA image", "FA image", us::Any(), true);
parser.addArgument("minFA", "fa", ctkCommandLineParser::Float, "Min. FA threshold", "minimum fractional anisotropy threshold", 0.15, true);
parser.addArgument("minCurv", "c", ctkCommandLineParser::Float, "Min. curvature radius", "minimum curvature radius in mm (default = 0.5*minimum-spacing)");
parser.addArgument("stepSize", "s", ctkCommandLineParser::Float, "Step size", "step size in mm (default = 0.1*minimum-spacing)");
parser.addArgument("tendf", "f", ctkCommandLineParser::Float, "Weight f", "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0, true);
parser.addArgument("tendg", "g", ctkCommandLineParser::Float, "Weight g", "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0, true);
parser.addArgument("numSeeds", "n", ctkCommandLineParser::Int, "Seeds per voxel", "Number of seeds per voxel.", 1, true);
parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "Min. fiber length", "minimum fiber length in mm", 20, true);
parser.addArgument("interpolate", "ip", ctkCommandLineParser::Bool, "Interpolate", "Use linear interpolation", false, true);
parser.addArgument("outFile", "o", ctkCommandLineParser::String, "Output file", "output fiber bundle (.fib)", us::Any(), false);
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setTitle("Streamline Tracking");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
ctkCommandLineParser::StringContainerType inputImages = us::any_cast<ctkCommandLineParser::StringContainerType>(parsedArgs["input"]);
string dtiFileName;
string outFileName = us::any_cast<string>(parsedArgs["outFile"]);
float minFA = 0.15;
float minCurv = -1;
float stepSize = -1;
float tendf = 1;
float tendg = 0;
float minLength = 20;
int numSeeds = 1;
bool interpolate = false;
if (parsedArgs.count("minCurv"))
minCurv = us::any_cast<float>(parsedArgs["minCurv"]);
if (parsedArgs.count("minFA"))
minFA = us::any_cast<float>(parsedArgs["minFA"]);
if (parsedArgs.count("stepSize"))
stepSize = us::any_cast<float>(parsedArgs["stepSize"]);
if (parsedArgs.count("tendf"))
tendf = us::any_cast<float>(parsedArgs["tendf"]);
if (parsedArgs.count("tendg"))
tendg = us::any_cast<float>(parsedArgs["tendg"]);
if (parsedArgs.count("minLength"))
minLength = us::any_cast<float>(parsedArgs["minLength"]);
if (parsedArgs.count("numSeeds"))
numSeeds = us::any_cast<int>(parsedArgs["numSeeds"]);
if (parsedArgs.count("interpolate"))
interpolate = us::any_cast<bool>(parsedArgs["interpolate"]);
try
{
typedef itk::StreamlineTrackingFilter< float > FilterType;
FilterType::Pointer filter = FilterType::New();
mitk::Image::Pointer mitkImage = NULL;
MITK_INFO << "Loading tensor images ...";
typedef itk::Image< itk::DiffusionTensor3D<float>, 3 > ItkTensorImage;
dtiFileName = inputImages.at(0);
for (unsigned int i=0; i<inputImages.size(); i++)
{
try
{
mitkImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(inputImages.at(i))->GetData());
mitk::TensorImage::Pointer img = dynamic_cast<mitk::TensorImage*>(mitk::IOUtil::LoadDataNode(inputImages.at(i))->GetData());
ItkTensorImage::Pointer itk_dti = ItkTensorImage::New();
mitk::CastToItkImage(img, itk_dti);
filter->SetInput(i, itk_dti);
}
catch(...){ MITK_INFO << "could not load: " << inputImages.at(i); }
}
MITK_INFO << "Loading seed image ...";
typedef itk::Image< unsigned char, 3 > ItkUCharImageType;
mitk::Image::Pointer mitkSeedImage = NULL;
if (parsedArgs.count("seed"))
mitkSeedImage = mitk::IOUtil::LoadImage(us::any_cast<string>(parsedArgs["seed"]));
MITK_INFO << "Loading mask image ...";
mitk::Image::Pointer mitkMaskImage = NULL;
if (parsedArgs.count("mask"))
mitkMaskImage = mitk::IOUtil::LoadImage(us::any_cast<string>(parsedArgs["mask"]));
// instantiate tracker
filter->SetSeedsPerVoxel(numSeeds);
filter->SetFaThreshold(minFA);
filter->SetMinCurvatureRadius(minCurv);
filter->SetStepSize(stepSize);
filter->SetF(tendf);
filter->SetG(tendg);
filter->SetInterpolate(interpolate);
filter->SetMinTractLength(minLength);
if (mitkSeedImage.IsNotNull())
{
ItkUCharImageType::Pointer mask = ItkUCharImageType::New();
mitk::CastToItkImage(mitkSeedImage, mask);
filter->SetSeedImage(mask);
}
if (mitkMaskImage.IsNotNull())
{
ItkUCharImageType::Pointer mask = ItkUCharImageType::New();
mitk::CastToItkImage(mitkMaskImage, mask);
filter->SetMaskImage(mask);
}
filter->Update();
vtkSmartPointer<vtkPolyData> fiberBundle = filter->GetFiberPolyData();
if ( fiberBundle->GetNumberOfLines()==0 )
{
MITK_INFO << "No fibers reconstructed. Check parametrization.";
return EXIT_FAILURE;
}
mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle);
fib->SetReferenceImage(mitkImage);
mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters();
for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it)
{
if ( (*it)->CanWriteBaseDataType(fib.GetPointer()) ) {
(*it)->SetFileName( outFileName.c_str() );
(*it)->DoWrite( fib.GetPointer() );
}
}
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
- MITK_INFO << "DONE";
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(StreamlineTracking);
diff --git a/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp b/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp
index 29a5323c70..6f4b744a57 100644
--- a/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp
@@ -1,195 +1,196 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkIOUtil.h>
#include "mitkImage.h"
#include <mitkImageCast.h>
#include "mitkITKImageImport.h"
#include <iostream>
#include <fstream>
#include <mitkTensorImage.h>
#include <mitkDiffusionImage.h>
#include "itkTensorDerivedMeasurementsFilter.h"
#include "itkDiffusionTensor3DReconstructionImageFilter.h"
#include "ctkCommandLineParser.h"
#include <itkImageFileWriter.h>
#include <itkNrrdImageIO.h>
#include <itkDiffusionTensor3D.h>
typedef short DiffusionPixelType;
typedef double TTensorPixelType;
static void ExtractMapsAndSave(mitk::TensorImage::Pointer tensorImage, std::string filename, std::string postfix = "")
{
mitk::Image* image = dynamic_cast<mitk::Image*> (tensorImage.GetPointer());
typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType;
typedef itk::Image< TensorPixelType, 3 > TensorImageType;
TensorImageType::Pointer itkvol = TensorImageType::New();
mitk::CastToItkImage(image, itkvol);
typedef itk::TensorDerivedMeasurementsFilter<TTensorPixelType> MeasurementsType;
MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New();
measurementsCalculator->SetInput(itkvol.GetPointer() );
mitk::Image::Pointer map = mitk::Image::New();
// FA
measurementsCalculator->SetMeasure(MeasurementsType::FA);
measurementsCalculator->Update();
map->InitializeByItk( measurementsCalculator->GetOutput() );
map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() );
mitk::IOUtil::SaveImage(map, filename + "_dti_FA" + postfix + ".nrrd");
// MD
measurementsCalculator->SetMeasure(MeasurementsType::MD);
measurementsCalculator->Update();
map->InitializeByItk( measurementsCalculator->GetOutput() );
map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() );
mitk::IOUtil::SaveImage(map, filename + "_dti_MD" + postfix + ".nrrd");
// AD
measurementsCalculator->SetMeasure(MeasurementsType::AD);
measurementsCalculator->Update();
map->InitializeByItk( measurementsCalculator->GetOutput() );
map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() );
mitk::IOUtil::SaveImage(map, filename + "_dti_AD" + postfix + ".nrrd");
// CA
measurementsCalculator->SetMeasure(MeasurementsType::CA);
measurementsCalculator->Update();
map->InitializeByItk( measurementsCalculator->GetOutput() );
map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() );
mitk::IOUtil::SaveImage(map, filename + "_dti_CA" + postfix + ".nrrd");
// RA
measurementsCalculator->SetMeasure(MeasurementsType::RA);
measurementsCalculator->Update();
map->InitializeByItk( measurementsCalculator->GetOutput() );
map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() );
mitk::IOUtil::SaveImage(map, filename + "_dti_RA" + postfix + ".nrrd");
// RD
measurementsCalculator->SetMeasure(MeasurementsType::RD);
measurementsCalculator->Update();
map->InitializeByItk( measurementsCalculator->GetOutput() );
map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() );
mitk::IOUtil::SaveImage(map, filename + "_dti_RD" + postfix + ".nrrd");
}
int TensorDerivedMapsExtraction(int argc, char* argv[])
{
+ MITK_INFO << "TensorDerivedMapsExtraction";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("help", "h", ctkCommandLineParser::String, "Help", "Show this help text");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input file", "input dwi file", us::Any(),false);
parser.addArgument("out", "o", ctkCommandLineParser::String, "Output folder", "output folder and base name, e.g. /tmp/outPatient1 ", us::Any(),false);
parser.setCategory("Diffusion Related Measures");
parser.setTitle("Tensor Derived Maps Extraction");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h"))
{
std::cout << "\n\nMiniApp Description: \nPerforms tensor reconstruction on DWI file," << endl;
std::cout << "and computes tensor derived measures." << endl;
std::cout << "\n\n For out parameter /tmp/outPatient1 it will produce :"<< endl;
std::cout << " /tmp/outPatient1_dti.dti , /tmp/outPatient1_dti_FA.nrrd, ..."<< endl;
std::cout << "\n\n Parameters:"<< endl;
std::cout << parser.helpText();
return EXIT_SUCCESS;
}
std::string inputFile = us::any_cast<string>(parsedArgs["input"]);
std::string baseFileName = us::any_cast<string>(parsedArgs["out"]);
std::string dtiFileName = "_dti.dti";
MITK_INFO << "BaseFileName: " << baseFileName;
mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputFile);
mitk::DiffusionImage<short>* diffusionImage = static_cast<mitk::DiffusionImage<short>*>(inputImage.GetPointer());
if (diffusionImage == NULL) // does NULL pointer check make sense after static cast ?
{
MITK_ERROR << "Invalid Input Image. Must be DWI. Aborting.";
return -1;
}
mitk::DiffusionImage<DiffusionPixelType>* vols = dynamic_cast <mitk::DiffusionImage<DiffusionPixelType>*> (inputImage.GetPointer());
typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType;
TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New();
typedef mitk::DiffusionImage<DiffusionPixelType> DiffusionImageType;
typedef DiffusionImageType::GradientDirectionContainerType GradientDirectionContainerType;
GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New();
for(GradientDirectionContainerType::ConstIterator it = vols->GetDirections()->Begin(); it != vols->GetDirections()->End(); it++)
{
gradientContainerCopy->push_back(it.Value());
}
tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() );
tensorReconstructionFilter->SetBValue(vols->GetReferenceBValue());
tensorReconstructionFilter->SetThreshold(50);
tensorReconstructionFilter->Update();
typedef itk::Image<itk::DiffusionTensor3D<TTensorPixelType>, 3> TensorImageType;
TensorImageType::Pointer tensorImage = tensorReconstructionFilter->GetOutput();
tensorImage->SetDirection( vols->GetVectorImage()->GetDirection() );
mitk::TensorImage::Pointer tensorImageMitk = mitk::TensorImage::New();
tensorImageMitk->InitializeByItk(tensorImage.GetPointer());
tensorImageMitk->SetVolume( tensorImage->GetBufferPointer() );
itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New();
io->SetFileType( itk::ImageIOBase::Binary );
io->UseCompressionOn();
itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< double >, 3 > >::Pointer writer = itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< double >, 3 > >::New();
writer->SetInput(tensorReconstructionFilter->GetOutput());
writer->SetFileName(baseFileName + dtiFileName);
writer->SetImageIO(io);
writer->UseCompressionOn();
writer->Update();
ExtractMapsAndSave(tensorImageMitk,baseFileName);
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(TensorDerivedMapsExtraction);
diff --git a/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp b/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp
index 446bcb6dfa..e73a4a4d01 100644
--- a/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp
@@ -1,102 +1,100 @@
/*===================================================================
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 "MiniAppManager.h"
#include "mitkBaseDataIOFactory.h"
#include "mitkDiffusionImage.h"
#include "mitkBaseData.h"
#include <itkDiffusionTensor3DReconstructionImageFilter.h>
#include <itkDiffusionTensor3D.h>
#include <itkImageFileWriter.h>
#include <itkNrrdImageIO.h>
#include "ctkCommandLineParser.h"
#include <itksys/SystemTools.hxx>
using namespace mitk;
/**
* Convert files from one ending to the other
*/
int TensorReconstruction(int argc, char* argv[])
{
+ MITK_INFO << "TensorReconstruction";
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::InputFile, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false);
parser.addArgument("outFile", "o", ctkCommandLineParser::OutputFile, "Output file", "output file", us::Any(), false);
parser.addArgument("b0Threshold", "t", ctkCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0, true);
parser.setCategory("Preprocessing Tools");
parser.setTitle("Tensor Reconstruction");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
std::string inFileName = us::any_cast<string>(parsedArgs["input"]);
std::string outfilename = us::any_cast<string>(parsedArgs["outFile"]);
outfilename = itksys::SystemTools::GetFilenamePath(outfilename)+"/"+itksys::SystemTools::GetFilenameWithoutExtension(outfilename);
outfilename += ".dti";
int threshold = 0;
if (parsedArgs.count("b0Threshold"))
threshold = us::any_cast<int>(parsedArgs["b0Threshold"]);
try
{
- MITK_INFO << "Loading image ...";
const std::string s1="", s2="";
std::vector<BaseData::Pointer> infile = BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false );
DiffusionImage<short>::Pointer dwi = dynamic_cast<DiffusionImage<short>*>(infile.at(0).GetPointer());
- MITK_INFO << "B0 threshold: " << threshold;
typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, float > TensorReconstructionImageFilterType;
TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New();
filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
filter->SetBValue(dwi->GetReferenceBValue());
filter->SetThreshold(threshold);
filter->Update();
// Save tensor image
- MITK_INFO << "writing image " << outfilename;
itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New();
io->SetFileType( itk::ImageIOBase::Binary );
io->UseCompressionOn();
itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< float >, 3 > >::Pointer writer = itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< float >, 3 > >::New();
writer->SetInput(filter->GetOutput());
writer->SetFileName(outfilename);
writer->SetImageIO(io);
writer->UseCompressionOn();
writer->Update();
}
catch ( itk::ExceptionObject &err)
{
MITK_INFO << "Exception: " << err;
}
catch ( std::exception err)
{
MITK_INFO << "Exception: " << err.what();
}
catch ( ... )
{
MITK_INFO << "Exception!";
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(TensorReconstruction);
diff --git a/Modules/DiffusionImaging/MiniApps/TractogramAngularError.cpp b/Modules/DiffusionImaging/MiniApps/TractogramAngularError.cpp
index 2b745330de..aa2a25ef82 100755
--- a/Modules/DiffusionImaging/MiniApps/TractogramAngularError.cpp
+++ b/Modules/DiffusionImaging/MiniApps/TractogramAngularError.cpp
@@ -1,205 +1,202 @@
/*===================================================================
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 "MiniAppManager.h"
#include <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <mitkDiffusionCoreObjectFactory.h>
#include <mitkFiberTrackingObjectFactory.h>
#include <metaCommand.h>
#include "ctkCommandLineParser.h"
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
#include <itkEvaluateTractogramDirectionsFilter.h>
#define _USE_MATH_DEFINES
#include <math.h>
int TractogramAngularError(int argc, char* argv[])
{
ctkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", us::Any(), false);
parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", us::Any(), false);
parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image");
parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results");
parser.addArgument("ignore", "n", ctkCommandLineParser::Bool, "don't increase error for missing or too many directions");
parser.addArgument("trilinear", "t", ctkCommandLineParser::Bool, "use trilinear instead of nearest neighbor interpolation");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
ctkCommandLineParser::StringContainerType referenceImages = us::any_cast<ctkCommandLineParser::StringContainerType>(parsedArgs["reference"]);
string fibFile = us::any_cast<string>(parsedArgs["input"]);
string maskImage("");
if (parsedArgs.count("mask"))
maskImage = us::any_cast<string>(parsedArgs["mask"]);
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;
if (parsedArgs.count("ignore"))
ignore = us::any_cast<bool>(parsedArgs["ignore"]);
bool interpolate = false;
if (parsedArgs.count("interpolate"))
interpolate = us::any_cast<bool>(parsedArgs["interpolate"]);
try
{
RegisterDiffusionCoreObjectFactory();
RegisterFiberTrackingObjectFactory();
typedef itk::Image<unsigned char, 3> ItkUcharImgType;
typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType;
typedef itk::VectorContainer< int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType;
typedef itk::EvaluateTractogramDirectionsFilter< float > EvaluationFilterType;
// load fiber bundle
mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast<mitk::FiberBundleX*>(mitk::IOUtil::LoadDataNode(fibFile)->GetData());
if (!inputTractogram)
return EXIT_FAILURE;
// load reference directions
ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New();
for (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(...){ MITK_INFO << "could not load: " << referenceImages.at(i); }
}
// load/create mask image
ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New();
if (maskImage.compare("")==0)
{
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);
}
else
{
mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(maskImage)->GetData());
mitk::CastToItkImage<ItkUcharImgType>(mitkMaskImage, itkMaskImage);
}
// evaluate directions
EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New();
evaluationFilter->SetTractogram(inputTractogram);
evaluationFilter->SetReferenceImageSet(referenceImageContainer);
evaluationFilter->SetMaskImage(itkMaskImage);
evaluationFilter->SetIgnoreMissingDirections(ignore);
evaluationFilter->SetUseInterpolation(interpolate);
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");
- MITK_INFO << "writing " << outfilename;
writer->SetFileName(outfilename.c_str());
writer->SetInput(angularErrorImage);
writer->Update();
}
string logFile = outRoot;
logFile.append("_ANGULAR_ERROR.csv");
ofstream file;
file.open (logFile.c_str());
string sens = "Mean:";
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMeanAngularError()));
sens.append(";\n");
sens.append("Median:");
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMedianAngularError()));
sens.append(";\n");
sens.append("Maximum:");
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMaxAngularError()));
sens.append(";\n");
sens.append("Minimum:");
sens.append(",");
sens.append(boost::lexical_cast<string>(evaluationFilter->GetMinAngularError()));
sens.append(";\n");
sens.append("STDEV:");
sens.append(",");
sens.append(boost::lexical_cast<string>(std::sqrt(evaluationFilter->GetVarAngularError())));
sens.append(";\n");
file << sens;
file.close();
-
- MITK_INFO << "DONE";
}
catch (itk::ExceptionObject e)
{
MITK_INFO << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
MITK_INFO << e.what();
return EXIT_FAILURE;
}
catch (...)
{
MITK_INFO << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
RegisterDiffusionMiniApp(TractogramAngularError);
diff --git a/Modules/DiffusionImaging/MiniApps/files.cmake b/Modules/DiffusionImaging/MiniApps/files.cmake
index 98428da34c..a2bcaa1eba 100644
--- a/Modules/DiffusionImaging/MiniApps/files.cmake
+++ b/Modules/DiffusionImaging/MiniApps/files.cmake
@@ -1,31 +1,32 @@
set(CPP_FILES
mitkDiffusionMiniApps.cpp
MiniAppManager.cpp
BatchedFolderRegistration.cpp
DicomFolderDump.cpp
FileFormatConverter.cpp
TensorReconstruction.cpp
TensorDerivedMapsExtraction.cpp
QballReconstruction.cpp
DiffusionIndices.cpp
ExtractImageStatistics.cpp
CopyGeometry.cpp
GibbsTracking.cpp
StreamlineTracking.cpp
FiberProcessing.cpp
LocalDirectionalFiberPlausibility.cpp
#TractogramAngularError.cpp
FiberDirectionExtraction.cpp
ImageResampler.cpp
PeakExtraction.cpp
PeaksAngularError.cpp
MultishellMethods.cpp
Fiberfox.cpp
ExportShImage.cpp
+ ImportShImage.cpp
NetworkCreation.cpp
NetworkStatistics.cpp
DwiDenoising.cpp
FiberExtraction.cpp
FiberJoin.cpp
)