diff --git a/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp b/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp
index de9023fb32..7fb1ae16ac 100755
--- a/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp
+++ b/Modules/DiffusionImaging/MiniApps/CopyGeometry.cpp
@@ -1,86 +1,81 @@
/*===================================================================
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 <mitkImageCast.h>
-#include <mitkDiffusionImage.h>
+#include <mitkImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkIOUtil.h>
#include "mitkCommandLineParser.h"
using namespace mitk;
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setTitle("Copy Geometry");
parser.setCategory("Preprocessing Tools");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("in", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(), false);
parser.addArgument("ref", "r", mitkCommandLineParser::InputFile, "Reference:", "reference image", us::Any(), false);
parser.addArgument("out", "o", mitkCommandLineParser::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
{
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::IOUtil::Save(dynamic_cast<DiffusionImage<short>*>(target.GetPointer()), outImage.c_str());
- }
- else
- mitk::IOUtil::SaveImage(target, outImage);
+ mitk::IOUtil::Save(target.GetPointer(), outImage.c_str());
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/DICOMLoader.cpp b/Modules/DiffusionImaging/MiniApps/DICOMLoader.cpp
index f62f43f2c2..e84237b798 100644
--- a/Modules/DiffusionImaging/MiniApps/DICOMLoader.cpp
+++ b/Modules/DiffusionImaging/MiniApps/DICOMLoader.cpp
@@ -1,272 +1,277 @@
/*===================================================================
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 "mitkBaseDataIOFactory.h"
-#include "mitkDiffusionImage.h"
+#include "mitkImage.h"
#include "mitkBaseData.h"
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkImageCast.h>
+#include <mitkITKImageImport.h>
#include <itkImageFileWriter.h>
#include <itkNrrdImageIO.h>
#include "mitkCommandLineParser.h"
#include <itksys/SystemTools.hxx>
#include <itksys/Directory.hxx>
#include "mitkDiffusionDICOMFileReader.h"
#include "mitkDICOMTagBasedSorter.h"
#include "mitkDICOMSortByTag.h"
#include "itkMergeDiffusionImagesFilter.h"
#include <mitkIOUtil.h>
static mitk::StringList& GetInputFilenames()
{
static mitk::StringList inputs;
return inputs;
}
void SetInputFileNames( std::string input_directory )
{
// I. Get all files in directory
itksys::Directory input;
input.Load( input_directory.c_str() );
// II. Push back files
mitk::StringList inputlist;//, mergedlist;
for( unsigned long idx=0; idx<input.GetNumberOfFiles(); idx++)
{
if( ! itksys::SystemTools::FileIsDirectory( input.GetFile(idx)) )
{
std::string fullpath = input_directory + "/" + std::string( input.GetFile(idx) );
inputlist.push_back( itksys::SystemTools::ConvertToOutputPath( fullpath.c_str() ) );
MITK_INFO("dicom.loader.inputdir.addfile") << input.GetFile(idx);
}
}
/*mergedlist.reserve( GetInputFilenames().size() + inputlist.size() );
mergedlist.insert( mergedlist.end(), GetInputFilenames().begin(), GetInputFilenames().end() );
mergedlist.insert( mergedlist.end(), inputlist.begin(), inputlist.end() );*/
GetInputFilenames() = inputlist;//mergedlist;
MITK_INFO("dicom.loader.setinputfiles.end") << "[]";
}
-mitk::DiffusionImage<short>::Pointer ReadInDICOMFiles( mitk::StringList& input_files, std::string output_file )
+mitk::Image::Pointer ReadInDICOMFiles( mitk::StringList& input_files, std::string output_file )
{
// repeat test with some more realistic sorting
mitk::DiffusionDICOMFileReader::Pointer gdcmReader = mitk::DiffusionDICOMFileReader::New(); // this also tests destruction
mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New();
// Use tags as in Qmitk
// all the things that split by tag in DicomSeriesReader
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0028, 0x0010) ); // Number of Rows
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0028, 0x0011) ); // Number of Columns
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0028, 0x0030) ); // Pixel Spacing
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0018, 0x1164) ); // Imager Pixel Spacing
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0020, 0x0037) ); // Image Orientation (Patient) // TODO add tolerance parameter (l. 1572 of original code)
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0018, 0x0050) ); // Slice Thickness
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0028, 0x0008) ); // Number of Frames
tagSorter->AddDistinguishingTag( mitk::DICOMTag(0x0020, 0x0052) ); // Frame of Reference UID
mitk::DICOMSortCriterion::ConstPointer sorting =
mitk::DICOMSortByTag::New( mitk::DICOMTag(0x0020, 0x0013), // instance number
mitk::DICOMSortByTag::New( mitk::DICOMTag(0x0020, 0x0012) //acquisition number
).GetPointer()
).GetPointer();
tagSorter->SetSortCriterion( sorting );
MITK_INFO("dicom.loader.read.init") << "[]" ;
MITK_INFO("dicom.loader.read.inputs") << " " << input_files.size();
gdcmReader->SetInputFiles( input_files );
gdcmReader->AddSortingElement( tagSorter );
gdcmReader->AnalyzeInputFiles();
gdcmReader->LoadImages();
mitk::Image::Pointer loaded_image = gdcmReader->GetOutput(0).GetMitkImage();
- mitk::DiffusionImage<short>::Pointer d_img = static_cast<mitk::DiffusionImage<short>*>( loaded_image.GetPointer() );
-
- return d_img;
+ return loaded_image;
}
typedef short DiffusionPixelType;
typedef itk::VectorImage<DiffusionPixelType,3> DwiImageType;
typedef DwiImageType::PixelType DwiPixelType;
typedef DwiImageType::RegionType DwiRegionType;
typedef std::vector< DwiImageType::Pointer > DwiImageContainerType;
-typedef mitk::DiffusionImage<DiffusionPixelType> DiffusionImageType;
-typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType;
+typedef mitk::Image DiffusionImageType;
+typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientContainerType;
typedef std::vector< GradientContainerType::Pointer > GradientListContainerType;
void SearchForInputInSubdirs( std::string root_directory, std::string subdir_prefix , std::vector<DiffusionImageType::Pointer>& output_container)
{
// I. Get all dirs in directory
itksys::Directory rootdir;
rootdir.Load( root_directory.c_str() );
MITK_INFO("dicom.loader.setinputdirs.start") << "Prefix = " << subdir_prefix;
for( unsigned int idx=0; idx<rootdir.GetNumberOfFiles(); idx++)
{
std::string current_path = rootdir.GetFile(idx);
std::string directory_path = std::string(rootdir.GetPath()) + std::string("/") + current_path;
MITK_INFO("dicom.loader.inputrootdir.test") << "ProbePath: " << current_path;
MITK_INFO("dicom.loader.inputrootdir.test") << "IsDirectory: " << itksys::SystemTools::FileIsDirectory( itksys::SystemTools::ConvertToOutputPath( directory_path.c_str() ).c_str() )
<< " StartsWith: " << itksys::SystemTools::StringStartsWith( current_path.c_str(), subdir_prefix.c_str() );
// test for prefix
if( itksys::SystemTools::FileIsDirectory( itksys::SystemTools::ConvertToOutputPath( directory_path.c_str() ).c_str() )
&& itksys::SystemTools::StringStartsWith( current_path.c_str(), subdir_prefix.c_str() )
)
{
MITK_INFO("dicom.loader.inputrootdir.searchin") << directory_path;
SetInputFileNames( itksys::SystemTools::ConvertToOutputPath( directory_path.c_str() ) );
MITK_INFO("dicom.loader.inputrootdir.preload") << "[]" ;
DiffusionImageType::Pointer dwi = ReadInDICOMFiles( GetInputFilenames(), "" );
output_container.push_back( dwi );
}
}
MITK_INFO("dicom.loader.setinputdirs.end") << "[]";
}
using namespace mitk;
/**
* Read DICOM Files through the new (refactored) Diffusion DICOM Loader. It serves also as a first test of the new functionality, it will replace the
* current loading mechanism after the necessary parts are merged into the master branch.
*/
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("inputdir", "i", mitkCommandLineParser::String, "Input Directory" ,"input directory containing dicom files", us::Any(), false);
parser.addArgument("output", "o", mitkCommandLineParser::String, "Output File Name", "output file", us::Any(), false);
parser.addArgument("dwprefix", "p", mitkCommandLineParser::String, "Recursive Scan Prefix", "prefix for subfolders search rootdir is specified by the 'inputdir' argument value", us::Any(), true);
parser.addArgument("dryrun", "-s", mitkCommandLineParser::Bool, "Dry run","do not read, only look for input files ", us::Any(), true );
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
{
return EXIT_FAILURE;
}
std::string inputDirectory = us::any_cast<std::string>( parsedArgs["inputdir"] );
MITK_INFO << "Loading data from directory: " << inputDirectory;
// retrieve the prefix flag (if set)
bool search_for_subdirs = false;
std::string subdir_prefix;
if( parsedArgs.count("dwprefix"))
{
MITK_INFO << "Prefix specified, will search for subdirs in the input directory!";
subdir_prefix = us::any_cast<std::string>( parsedArgs["dwprefix"] );
search_for_subdirs = true;
}
// retrieve the output
std::string outputFile = us::any_cast< std::string >( parsedArgs["output"] );
// if the executable is called with a single directory, just parse the given folder for files and read them into a diffusion image
if( !search_for_subdirs )
{
SetInputFileNames( inputDirectory );
MITK_INFO << "Got " << GetInputFilenames().size() << " input files.";
- mitk::DiffusionImage<short>::Pointer d_img = ReadInDICOMFiles( GetInputFilenames(), outputFile );
+ mitk::Image::Pointer d_img = ReadInDICOMFiles( GetInputFilenames(), outputFile );
try
{
mitk::IOUtil::Save(d_img, outputFile.c_str());
}
catch( const itk::ExceptionObject& e)
{
MITK_ERROR << "Failed to write out the output file. \n\t Reason : ITK Exception " << e.what();
}
}
// if the --dwprefix flag is set, then we have to look for the directories, load each of them separately and afterwards merge the images
else
{
- std::vector<mitk::DiffusionImage<DiffusionPixelType>::Pointer> output_container;
+ std::vector<mitk::Image::Pointer> output_container;
SearchForInputInSubdirs( inputDirectory, subdir_prefix, output_container );
// final output image
- mitk::DiffusionImage<DiffusionPixelType>::Pointer image = mitk::DiffusionImage<DiffusionPixelType>::New();
+ mitk::Image::Pointer image = mitk::Image::New();
if( output_container.size() > 1 )
{
DwiImageContainerType imageContainer;
GradientListContainerType gradientListContainer;
std::vector< double > bValueContainer;
- for ( std::vector< mitk::DiffusionImage<DiffusionPixelType>::Pointer >::iterator dwi = output_container.begin();
- dwi != output_container.end(); ++dwi )
+ for ( std::vector< mitk::Image::Pointer >::iterator dwi = output_container.begin();
+ dwi != output_container.end(); ++dwi )
{
- imageContainer.push_back((*dwi)->GetVectorImage());
- gradientListContainer.push_back((*dwi)->GetDirections());
- bValueContainer.push_back((*dwi)->GetReferenceBValue());
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(*dwi, itkVectorImagePointer);
+
+ imageContainer.push_back(itkVectorImagePointer);
+ gradientListContainer.push_back( mitk::DiffusionPropertyHelper::GetGradientContainer(*dwi));
+ bValueContainer.push_back( mitk::DiffusionPropertyHelper::GetReferenceBValue(*dwi));
}
typedef itk::MergeDiffusionImagesFilter<short> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetImageVolumes(imageContainer);
filter->SetGradientLists(gradientListContainer);
filter->SetBValues(bValueContainer);
filter->Update();
vnl_matrix_fixed< double, 3, 3 > mf; mf.set_identity();
- image->SetVectorImage( filter->GetOutput() );
- image->SetReferenceBValue(filter->GetB_Value());
- image->SetDirections(filter->GetOutputGradients());
- image->SetMeasurementFrame(mf);
- image->InitializeFromVectorImage();
+ image = mitk::GrabItkImageMemory( filter->GetOutput() );
+ image->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( filter->GetOutputGradients() ) );
+ image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( filter->GetB_Value() ) );
+ image->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( mf ) );
+ mitk::DiffusionPropertyHelper propertyHelper( image );
+ propertyHelper.InitializeImage();
}
// just output the image if there was only one folder found
else
{
image = output_container.at(0);
}
MITK_INFO("dicom.import.writeout") << " [OutputFile] " << outputFile.c_str();
try
{
mitk::IOUtil::Save(image, outputFile.c_str());
}
catch( const itk::ExceptionObject& e)
{
MITK_ERROR << "Failed to write out the output file. \n\t Reason : ITK Exception " << e.what();
}
}
return 1;
}
diff --git a/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp b/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp
index 83f8f18544..f30bfd5189 100644
--- a/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp
+++ b/Modules/DiffusionImaging/MiniApps/DwiDenoising.cpp
@@ -1,156 +1,163 @@
/*===================================================================
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 <mitkImageCast.h>
#include <mitkBaseDataIOFactory.h>
#include "mitkCommandLineParser.h"
#include <boost/algorithm/string.hpp>
-#include <DiffusionWeightedImages/mitkDiffusionImage.h>
+#include <mitkImage.h>
#include <itkNonLocalMeansDenoisingFilter.h>
#include <itkImage.h>
#include <mitkIOUtil.h>
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkImageCast.h>
+#include <mitkITKImageImport.h>
+#include <mitkProperties.h>
-typedef mitk::DiffusionImage<short> DiffusionImageType;
+typedef mitk::Image 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() )
{
std::cout << "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 main(int argc, char* argv[])
{
std::cout << "DwiDenoising";
mitkCommandLineParser 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", mitkCommandLineParser::InputFile, "Input:", "input image (DWI)", us::Any(), false);
parser.addArgument("variance", "v", mitkCommandLineParser::Float, "Variance:", "noise variance", us::Any(), false);
parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "brainmask for input image", us::Any(), true);
parser.addArgument("search", "s", mitkCommandLineParser::Int, "Search radius:", "search radius", us::Any(), true);
parser.addArgument("compare", "c", mitkCommandLineParser::Int, "Comparison radius:", "comparison radius", us::Any(), true);
parser.addArgument("joint", "j", mitkCommandLineParser::Bool, "Joint information:", "use joint information");
parser.addArgument("rician", "r", mitkCommandLineParser::Bool, "Rician adaption:", "use rician adaption");
parser.changeParameterGroup("Output", "Output of this miniapp");
parser.addArgument("output", "o", mitkCommandLineParser::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());
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
itk::NonLocalMeansDenoisingFilter<short>::Pointer filter = itk::NonLocalMeansDenoisingFilter<short>::New();
filter->SetNumberOfThreads(12);
- filter->SetInputImage(dwi->GetVectorImage());
+ filter->SetInputImage( itkVectorImagePointer );
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();
+ DiffusionImageType::Pointer output = mitk::GrabItkImageMemory( filter->GetOutput() );
+ output->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi) ) );
+ output->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi) ) );
+ mitk::DiffusionPropertyHelper propertyHelper( output );
+ propertyHelper.InitializeImage();
// std::stringstream name;
// name << outFileName << "_NLM_" << search << "-" << compare << "-" << variance << ".dwi";
mitk::IOUtil::Save(output, outFileName.c_str());
}
else
{
std::cout << "Only supported for .dwi!";
}
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp b/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
index 3bf5e834b8..c50a85159d 100755
--- a/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/FiberDirectionExtraction.cpp
@@ -1,174 +1,174 @@
/*===================================================================
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 <mitkBaseDataIOFactory.h>
#include <mitkBaseData.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <metaCommand.h>
#include "mitkCommandLineParser.h"
#include <usAny.h>
#include <itkImageFileWriter.h>
#include <mitkIOUtil.h>
#include <boost/lexical_cast.hpp>
#include <itkEvaluateDirectionImagesFilter.h>
#include <itkTractsToVectorImageFilter.h>
#include <mitkCoreObjectFactory.h>
#define _USE_MATH_DEFINES
#include <math.h>
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setTitle("Fiber Direction Extraction");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib/.trk)", us::Any(), false);
parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image");
parser.addArgument("athresh", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true);
parser.addArgument("peakthresh", "t", mitkCommandLineParser::Float, "Peak size threshold:", "peak size threshold relative to largest peak in voxel", 0.2, true);
parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results");
parser.addArgument("numdirs", "d", mitkCommandLineParser::Int, "Max. num. directions:", "maximum number of fibers per voxel", 3, true);
parser.addArgument("normalize", "n", mitkCommandLineParser::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)
{
std::cout << "Using mask image";
itkMaskImage = ItkUcharImgType::New();
mitk::Image::Pointer mitkMaskImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(maskImage)->GetData());
- mitk::CastToItkImage<ItkUcharImgType>(mitkMaskImage, itkMaskImage);
+ mitk::CastToItkImage(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");
writer->SetFileName(outfilename.c_str());
writer->SetInput(itkImg);
writer->Update();
}
if (verbose)
{
// write vector field
mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle();
string outfilename = outRoot;
outfilename.append("_VECTOR_FIELD.fib");
mitk::IOUtil::SaveBaseData(directions.GetPointer(), outfilename );
// write num direction image
{
ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage();
typedef itk::ImageFileWriter< ItkUcharImgType > WriterType;
WriterType::Pointer writer = WriterType::New();
string outfilename = outRoot;
outfilename.append("_NUM_DIRECTIONS.nrrd");
writer->SetFileName(outfilename.c_str());
writer->SetInput(numDirImage);
writer->Update();
}
}
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp b/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp
index e767f81004..324fecfba2 100755
--- a/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp
+++ b/Modules/DiffusionImaging/MiniApps/Fiberfox.cpp
@@ -1,78 +1,80 @@
/*===================================================================
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 <mitkImageCast.h>
-#include <mitkDiffusionImage.h>
+#include <mitkITKImageImport.h>
+#include <mitkProperties.h>
+#include <mitkImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkIOUtil.h>
#include <mitkFiberBundleX.h>
#include <mitkFiberfoxParameters.h>
#include "mitkCommandLineParser.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 */
using namespace mitk;
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output root:", "output root", us::Any(), false);
parser.addArgument("parameters", "p", mitkCommandLineParser::InputFile, "Parameter file:", "fiberfox parameter file", us::Any(), false);
parser.addArgument("fiberbundle", "f", mitkCommandLineParser::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();
+ mitk::Image::Pointer image = mitk::GrabItkImageMemory( tractsToDwiFilter->GetOutput() );
+ image->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() ) );
+ image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.m_Bvalue ) );
+ mitk::DiffusionPropertyHelper propertyHelper( image );
+ propertyHelper.InitializeImage();
mitk::IOUtil::Save(image, outName.c_str());
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp b/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp
index faff7c6965..a05ce2419b 100644
--- a/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp
+++ b/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp
@@ -1,326 +1,314 @@
/*===================================================================
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 "mitkCommandLineParser.h"
#include <mitkIOUtil.h>
#include <mitkRegistrationWrapper.h>
#include <mitkImage.h>
#include <mitkImageCast.h>
#include <mitkITKImageImport.h>
-#include <mitkDiffusionImage.h>
#include <mitkImageTimeSelector.h>
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkProperties.h>
// ITK
#include <itksys/SystemTools.hxx>
#include <itkDirectory.h>
#include "itkLinearInterpolateImageFunction.h"
#include "itkWindowedSincInterpolateImageFunction.h"
#include "itkIdentityTransform.h"
#include "itkResampleImageFilter.h"
#include "itkResampleDwiImageFilter.h"
typedef itk::Image<double, 3> InputImageType;
-typedef mitk::DiffusionImage<short> DiffusionImageType;
static mitk::Image::Pointer TransformToReference(mitk::Image *reference, mitk::Image *moving, bool sincInterpol = false)
{
// Convert to itk Images
InputImageType::Pointer itkReference = InputImageType::New();
InputImageType::Pointer itkMoving = InputImageType::New();
mitk::CastToItkImage(reference,itkReference);
mitk::CastToItkImage(moving,itkMoving);
// Identify Transform
typedef itk::IdentityTransform<double, 3> T_Transform;
T_Transform::Pointer _pTransform = T_Transform::New();
_pTransform->SetIdentity();
typedef itk::WindowedSincInterpolateImageFunction< InputImageType, 3> WindowedSincInterpolatorType;
WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New();
typedef itk::ResampleImageFilter<InputImageType, InputImageType> ResampleFilterType;
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetInput(itkMoving);
resampler->SetReferenceImage( itkReference );
resampler->UseReferenceImageOn();
resampler->SetTransform(_pTransform);
resampler->SetInterpolator(sinc_interpolator);
resampler->Update();
// Convert back to mitk
mitk::Image::Pointer result = mitk::Image::New();
result->InitializeByItk(resampler->GetOutput());
GrabItkImageMemory( resampler->GetOutput() , result );
return result;
}
static std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems)
{
std::stringstream ss(s);
std::string item;
while (std::getline(ss, item, delim))
{
elems.push_back(item);
}
return elems;
}
static std::vector<std::string> split(const std::string &s, char delim)
{
std::vector < std::string > elems;
return split(s, delim, elems);
}
static mitk::Image::Pointer ResampleBySpacing(mitk::Image *input, float *spacing, bool useLinInt = true)
{
InputImageType::Pointer itkImage = InputImageType::New();
CastToItkImage(input,itkImage);
/**
* 1) Resampling
*
*/
// Identity transform.
// We don't want any transform on our image except rescaling which is not
// specified by a transform but by the input/output spacing as we will see
// later.
// So no transform will be specified.
typedef itk::IdentityTransform<double, 3> T_Transform;
// The resampler type itself.
typedef itk::ResampleImageFilter<InputImageType, InputImageType> T_ResampleFilter;
// Prepare the resampler.
// Instantiate the transform and specify it should be the id transform.
T_Transform::Pointer _pTransform = T_Transform::New();
_pTransform->SetIdentity();
// Instantiate the resampler. Wire in the transform and the interpolator.
T_ResampleFilter::Pointer _pResizeFilter = T_ResampleFilter::New();
// Specify the input.
_pResizeFilter->SetInput(itkImage);
_pResizeFilter->SetTransform(_pTransform);
// Set the output origin.
_pResizeFilter->SetOutputOrigin(itkImage->GetOrigin());
// Compute the size of the output.
// The size (# of pixels) in the output is recomputed using
// the ratio of the input and output sizes.
InputImageType::SpacingType inputSpacing = itkImage->GetSpacing();
InputImageType::SpacingType outputSpacing;
const InputImageType::RegionType& inputSize = itkImage->GetLargestPossibleRegion();
InputImageType::SizeType outputSize;
typedef InputImageType::SizeType::SizeValueType SizeValueType;
// Set the output spacing.
outputSpacing[0] = spacing[0];
outputSpacing[1] = spacing[1];
outputSpacing[2] = spacing[2];
outputSize[0] = static_cast<SizeValueType>(inputSize.GetSize()[0] * inputSpacing[0] / outputSpacing[0] + .5);
outputSize[1] = static_cast<SizeValueType>(inputSize.GetSize()[1] * inputSpacing[1] / outputSpacing[1] + .5);
outputSize[2] = static_cast<SizeValueType>(inputSize.GetSize()[2] * inputSpacing[2] / outputSpacing[2] + .5);
_pResizeFilter->SetOutputSpacing(outputSpacing);
_pResizeFilter->SetSize(outputSize);
typedef itk::LinearInterpolateImageFunction< InputImageType > LinearInterpolatorType;
LinearInterpolatorType::Pointer lin_interpolator = LinearInterpolatorType::New();
typedef itk::WindowedSincInterpolateImageFunction< InputImageType, 4> WindowedSincInterpolatorType;
WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New();
if (useLinInt)
_pResizeFilter->SetInterpolator(lin_interpolator);
else
_pResizeFilter->SetInterpolator(sinc_interpolator);
_pResizeFilter->Update();
mitk::Image::Pointer image = mitk::Image::New();
image->InitializeByItk(_pResizeFilter->GetOutput());
mitk::GrabItkImageMemory( _pResizeFilter->GetOutput(), image);
return image;
}
/// Save images according to file type
static void SaveImage(std::string fileName, mitk::Image* image, std::string fileType )
{
std::cout << "----Save to " << fileName;
- if (fileType == "dwi") // IOUtil does not handle dwi files properly Bug 15772
- {
- try
- {
- mitk::IOUtil::Save(dynamic_cast<mitk::DiffusionImage<short>*>(image), fileName.c_str());
- }
- catch( const itk::ExceptionObject& e)
- {
- MITK_ERROR << "Caught exception: " << e.what();
- mitkThrow() << "Failed with exception from subprocess!";
- }
- }
- else
- {
- mitk::IOUtil::SaveImage(image, fileName);
- }
+ mitk::IOUtil::Save(image, fileName);
}
-DiffusionImageType::Pointer ResampleDWIbySpacing(DiffusionImageType::Pointer input, float* spacing, bool useLinInt = true)
+mitk::Image::Pointer ResampleDWIbySpacing(mitk::Image::Pointer input, float* spacing, bool useLinInt = true)
{
itk::Vector<double, 3> spacingVector;
spacingVector[0] = spacing[0];
spacingVector[1] = spacing[1];
spacingVector[2] = spacing[2];
typedef itk::ResampleDwiImageFilter<short> ResampleFilterType;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(input, itkVectorImagePointer);
+
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
- resampler->SetInput(input->GetVectorImage());
+ resampler->SetInput( itkVectorImagePointer );
resampler->SetInterpolation(ResampleFilterType::Interpolate_Linear);
resampler->SetNewSpacing(spacingVector);
resampler->Update();
- DiffusionImageType::Pointer output = DiffusionImageType::New();
- output->SetVectorImage(resampler->GetOutput());
- output->SetDirections(input->GetDirections());
- output->SetReferenceBValue(input->GetReferenceBValue());
- output->InitializeFromVectorImage();
+ mitk::Image::Pointer output = mitk::GrabItkImageMemory( resampler->GetOutput() );
+ output->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( mitk::DiffusionPropertyHelper::GetGradientContainer(input) ) );
+ output->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( mitk::DiffusionPropertyHelper::GetReferenceBValue(input) ) );
+ mitk::DiffusionPropertyHelper propertyHelper( output );
+ propertyHelper.InitializeImage();
return output;
}
int main( int argc, char* argv[] )
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--","-");
parser.setTitle("Image Resampler");
parser.setCategory("Preprocessing Tools");
parser.setContributor("MBI");
parser.setDescription("Resample an image to eigther a specific spacing or to a reference image.");
// Add command line argument names
parser.addArgument("help", "h",mitkCommandLineParser::Bool, "Show this help text");
parser.addArgument("input", "i", mitkCommandLineParser::InputImage, "Input:", "Input file",us::Any(),false);
parser.addArgument("output", "o", mitkCommandLineParser::OutputDirectory, "Output:", "Output folder (ending with /)",us::Any(),false);
parser.addArgument("spacing", "s", mitkCommandLineParser::String, "Spacing:", "Resample provide x,y,z spacing in mm (e.g. -r 1,1,3), is not applied to tensor data",us::Any());
parser.addArgument("reference", "r", mitkCommandLineParser::String, "Reference:", "Resample using supplied reference image. Also cuts image to same dimensions",us::Any());
parser.addArgument("win-sinc", "w", mitkCommandLineParser::Bool, "Windowed-sinc interpolation:", "Use windowed-sinc interpolation (3) instead of linear interpolation ",us::Any());
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
// Handle special arguments
bool useSpacing = false;
bool useLinearInterpol = true;
{
if (parsedArgs.size() == 0)
{
return EXIT_FAILURE;
}
if (parsedArgs.count("sinc-int"))
useLinearInterpol = false;
// Show a help message
if ( parsedArgs.count("help") || parsedArgs.count("h"))
{
std::cout << parser.helpText();
return EXIT_SUCCESS;
}
}
std::string outputPath = us::any_cast<string>(parsedArgs["output"]);
std::string inputFile = us::any_cast<string>(parsedArgs["input"]);
std::vector<std::string> spacings;
float spacing[3];
if (parsedArgs.count("spacing"))
{
std::string arg = us::any_cast<string>(parsedArgs["spacing"]);
spacings = split(arg ,',');
spacing[0] = atoi(spacings.at(0).c_str());
spacing[1] = atoi(spacings.at(1).c_str());
spacing[2] = atoi(spacings.at(2).c_str());
useSpacing = true;
}
std::string refImageFile = "";
if (parsedArgs.count("reference"))
{
refImageFile = us::any_cast<string>(parsedArgs["reference"]);
}
if (refImageFile =="" && useSpacing == false)
{
MITK_ERROR << "No information how to resample is supplied. Use eigther --spacing or --reference !";
return EXIT_FAILURE;
}
mitk::Image::Pointer refImage;
if (!useSpacing)
refImage = mitk::IOUtil::LoadImage(refImageFile);
- DiffusionImageType::Pointer inputDWI = dynamic_cast<DiffusionImageType*>(mitk::IOUtil::LoadBaseData(inputFile).GetPointer());
- if (inputDWI.IsNotNull())
+ mitk::Image::Pointer inputDWI = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadBaseData(inputFile).GetPointer());
+ if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(inputDWI))
{
- DiffusionImageType::Pointer outputImage;
+ mitk::Image::Pointer outputImage;
if (useSpacing)
outputImage = ResampleDWIbySpacing(inputDWI, spacing);
else
{
MITK_WARN << "Not supported yet, to resample a DWI please set a new spacing.";
return EXIT_FAILURE;
}
std::string fileStem = itksys::SystemTools::GetFilenameWithoutExtension(inputFile);
std::string outName(outputPath + fileStem + "_res.dwi");
mitk::IOUtil::Save(outputImage, outName.c_str());
return EXIT_SUCCESS;
}
mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputFile);
mitk::Image::Pointer resultImage;
if (useSpacing)
resultImage = ResampleBySpacing(inputImage,spacing);
else
resultImage = TransformToReference(refImage,inputImage);
std::string fileStem = itksys::SystemTools::GetFilenameWithoutExtension(inputFile);
mitk::IOUtil::SaveImage(resultImage, outputPath + fileStem + "_res.nrrd");
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp b/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp
index 687ed57bf1..bc81479de8 100644
--- a/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp
+++ b/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp
@@ -1,215 +1,219 @@
/*===================================================================
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 <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 <mitkImage.h>
#include <mitkQBallImage.h>
#include <mitkBaseData.h>
#include <mitkFiberBundleX.h>
#include "mitkCommandLineParser.h"
#include <boost/lexical_cast.hpp>
#include <itkRadialMultishellToSingleshellImageFilter.h>
#include <itkADCAverageFunctor.h>
#include <itkBiExpFitFunctor.h>
#include <itkKurtosisFitFunctor.h>
#include <itkDwiGradientLengthCorrectionFilter.h>
#include <mitkIOUtil.h>
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkProperties.h>
+#include <mitkImageCast.h>
+#include <mitkITKImageImport.h>
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setTitle("Multishell Methods");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("in", "i", mitkCommandLineParser::InputFile, "Input:", "input file", us::Any(), false);
parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output:", "output file", us::Any(), false);
parser.addArgument("adc", "D", mitkCommandLineParser::Bool, "ADC:", "ADC Average", us::Any(), false);
parser.addArgument("akc", "K", mitkCommandLineParser::Bool, "Kurtosis fit:", "Kurtosis Fit", us::Any(), false);
parser.addArgument("biexp", "B", mitkCommandLineParser::Bool, "BiExp fit:", "BiExp fit", us::Any(), false);
parser.addArgument("targetbvalue", "b", mitkCommandLineParser::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
{
std::cout << "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()) )
+ if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast<mitk::Image*>(baseData.GetPointer()) ) )
{
- mitk::DiffusionImage<short>::Pointer dwi = dynamic_cast<mitk::DiffusionImage<short>*>(baseData.GetPointer());
+ mitk::Image::Pointer dwi = dynamic_cast<mitk::Image*>(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->SetReferenceBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi ));
+ roundfilter->SetReferenceGradientDirectionContainer(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi));
roundfilter->Update();
- dwi->SetReferenceBValue( roundfilter->GetNewBValue() );
- dwi->SetDirections( roundfilter->GetOutputGradientDirectionContainer());
+ dwi->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( roundfilter->GetNewBValue() ) );
+ dwi->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( roundfilter->GetOutputGradientDirectionContainer() ) );
// filter input parameter
- const mitk::DiffusionImage<short>::BValueMap
- &originalShellMap = dwi->GetBValueMap();
+ const mitk::DiffusionPropertyHelper::BValueMapType
+ &originalShellMap = mitk::DiffusionPropertyHelper::GetBValueMap(dwi);
- const mitk::DiffusionImage<short>::ImageType
- *vectorImage = dwi->GetVectorImage();
+ mitk::DiffusionPropertyHelper::ImageType::Pointer vectorImage = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, vectorImage);
- const mitk::DiffusionImage<short>::GradientDirectionContainerType::Pointer
- gradientContainer = dwi->GetDirections();
+ const mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer
+ gradientContainer = mitk::DiffusionPropertyHelper::GetGradientContainer(dwi);
const unsigned int
- &bValue = dwi->GetReferenceBValue();
+ &bValue = mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi );
// filter call
vnl_vector<double> bValueList(originalShellMap.size()-1);
double targetBValue = bValueList.mean();
- mitk::DiffusionImage<short>::BValueMap::const_iterator it = originalShellMap.begin();
+ mitk::DiffusionPropertyHelper::BValueMapType::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::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( targetBValue ) );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( filter->GetTargetGradientDirections() ) );
+ mitk::DiffusionPropertyHelper propertyHelper( outImage );
+ propertyHelper.InitializeImage();
mitk::IOUtil::Save(outImage, (outName + "_ADC.dwi").c_str());
}
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::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( targetBValue ) );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( filter->GetTargetGradientDirections() ) );
+ mitk::DiffusionPropertyHelper propertyHelper( outImage );
+ propertyHelper.InitializeImage();
mitk::IOUtil::Save(outImage, (string(outName) + "_AKC.dwi").c_str());
}
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::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( targetBValue ) );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( filter->GetTargetGradientDirections() ) );
+ mitk::DiffusionPropertyHelper propertyHelper( outImage );
+ propertyHelper.InitializeImage();
mitk::IOUtil::Save(outImage, (string(outName) + "_BiExp.dwi").c_str());
}
}
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp
index 1fe492899c..2d124f8659 100644
--- a/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp
+++ b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp
@@ -1,515 +1,516 @@
/*===================================================================
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.
===================================================================*/
// std includes
#include <string>
#include <sstream>
#include <fstream>
#include <vector>
#include <map>
#include <utility>
// boost includes
#include <boost/algorithm/string.hpp>
// ITK includes
#include <itkImageFileWriter.h>
// CTK includes
#include "mitkCommandLineParser.h"
// MITK includes
#include <mitkBaseDataIOFactory.h>
#include <mitkConnectomicsStatisticsCalculator.h>
#include <mitkConnectomicsNetworkThresholder.h>
#include <itkConnectomicsNetworkToConnectivityMatrixImageFilter.h>
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
+
parser.addArgument("inputNetwork", "i", mitkCommandLineParser::InputFile, "Input network", "input connectomics network (.cnf)", us::Any(), false);
parser.addArgument("outputFile", "o", mitkCommandLineParser::OutputFile, "Output file", "name of output file", us::Any(), false);
parser.addArgument("noGlobalStatistics", "g", mitkCommandLineParser::Bool, "No global statistics", "Do not calculate global statistics");
parser.addArgument("createConnectivityMatriximage", "I", mitkCommandLineParser::Bool, "Write connectivity matrix image", "Write connectivity matrix image");
parser.addArgument("binaryConnectivity", "b", mitkCommandLineParser::Bool, "Binary connectivity", "Whether to create a binary connectivity matrix");
parser.addArgument("rescaleConnectivity", "r", mitkCommandLineParser::Bool, "Rescale connectivity", "Whether to rescale the connectivity matrix");
parser.addArgument("localStatistics", "L", mitkCommandLineParser::StringList, "Local statistics", "Provide a list of node labels for local statistics", us::Any());
parser.addArgument("regionList", "R", mitkCommandLineParser::StringList, "Region list", "A space separated list of regions. Each region has the format\n regionname;label1;label2;...;labelN", us::Any());
parser.addArgument("granularity", "gr", mitkCommandLineParser::Int, "Granularity", "How finely to test the density range and how many thresholds to consider");
parser.addArgument("startDensity", "d", mitkCommandLineParser::Float, "Start Density", "Largest density for the range");
parser.addArgument("thresholdStepSize", "t", mitkCommandLineParser::Int, "Step size threshold", "Distance of two adjacent thresholds");
parser.setCategory("Connectomics");
parser.setTitle("Network Statistics");
parser.setDescription("");
parser.setContributor("MBI");
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
//default values
bool noGlobalStatistics( false );
bool binaryConnectivity( false );
bool rescaleConnectivity( false );
bool createConnectivityMatriximage( false );
int granularity( 1 );
double startDensity( 1.0 );
int thresholdStepSize( 3 );
// parse command line arguments
std::string networkName = us::any_cast<std::string>(parsedArgs["inputNetwork"]);
std::string outName = us::any_cast<std::string>(parsedArgs["outputFile"]);
mitkCommandLineParser::StringContainerType localLabels;
if(parsedArgs.count("localStatistics"))
{
localLabels = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["localStatistics"]);
}
mitkCommandLineParser::StringContainerType unparsedRegions;
std::map< std::string, std::vector<std::string> > parsedRegions;
std::map< std::string, std::vector<std::string> >::iterator parsedRegionsIterator;
if(parsedArgs.count("regionList"))
{
unparsedRegions = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["regionList"]);
for(unsigned int index(0); index < unparsedRegions.size(); index++ )
{
std::vector< std::string > tempRegionVector;
boost::split(tempRegionVector, unparsedRegions.at(index), boost::is_any_of(";"));
std::vector< std::string >::const_iterator begin = tempRegionVector.begin();
std::vector< std::string >::const_iterator last = tempRegionVector.begin() + tempRegionVector.size();
std::vector< std::string > insertRegionVector(begin + 1, last);
if( parsedRegions.count( tempRegionVector.at(0) ) == 0 )
{
parsedRegions.insert( std::pair< std::string, std::vector<std::string> >( tempRegionVector.at(0), insertRegionVector) );
}
else
{
MITK_ERROR << "Region already exists. Skipping second occurrence.";
}
}
}
if (parsedArgs.count("noGlobalStatistics"))
noGlobalStatistics = us::any_cast<bool>(parsedArgs["noGlobalStatistics"]);
if (parsedArgs.count("binaryConnectivity"))
binaryConnectivity = us::any_cast<bool>(parsedArgs["binaryConnectivity"]);
if (parsedArgs.count("rescaleConnectivity"))
rescaleConnectivity = us::any_cast<bool>(parsedArgs["rescaleConnectivity"]);
if (parsedArgs.count("createConnectivityMatriximage"))
createConnectivityMatriximage = us::any_cast<bool>(parsedArgs["createConnectivityMatriximage"]);
if (parsedArgs.count("granularity"))
granularity = us::any_cast<int>(parsedArgs["granularity"]);
if (parsedArgs.count("startDensity"))
startDensity = us::any_cast<float>(parsedArgs["startDensity"]);
if (parsedArgs.count("thresholdStepSize"))
thresholdStepSize = us::any_cast<int>(parsedArgs["thresholdStepSize"]);
try
{
const std::string s1="", s2="";
// load network
std::vector<mitk::BaseData::Pointer> networkFile =
mitk::BaseDataIO::LoadBaseDataFromFile( networkName, s1, s2, false );
if( networkFile.empty() )
{
std::string errorMessage = "File at " + networkName + " could not be read. Aborting.";
MITK_ERROR << errorMessage;
return EXIT_FAILURE;
}
mitk::BaseData* networkBaseData = networkFile.at(0);
mitk::ConnectomicsNetwork* network = dynamic_cast<mitk::ConnectomicsNetwork*>( networkBaseData );
if( !network )
{
std::string errorMessage = "Read file at " + networkName + " could not be recognized as network. Aborting.";
MITK_ERROR << errorMessage;
return EXIT_FAILURE;
}
// streams
std::stringstream globalHeaderStream;
globalHeaderStream << "NumberOfVertices "
<< "NumberOfEdges "
<< "AverageDegree "
<< "ConnectionDensity "
<< "NumberOfConnectedComponents "
<< "AverageComponentSize "
<< "LargestComponentSize "
<< "RatioOfNodesInLargestComponent "
<< "HopPlotExponent "
<< "EffectiveHopDiameter "
<< "AverageClusteringCoefficientsC "
<< "AverageClusteringCoefficientsD "
<< "AverageClusteringCoefficientsE "
<< "AverageVertexBetweennessCentrality "
<< "AverageEdgeBetweennessCentrality "
<< "NumberOfIsolatedPoints "
<< "RatioOfIsolatedPoints "
<< "NumberOfEndPoints "
<< "RatioOfEndPoints "
<< "Diameter "
<< "Diameter90 "
<< "Radius "
<< "Radius90 "
<< "AverageEccentricity "
<< "AverageEccentricity90 "
<< "AveragePathLength "
<< "NumberOfCentralPoints "
<< "RatioOfCentralPoints "
<< "SpectralRadius "
<< "SecondLargestEigenValue "
<< "AdjacencyTrace "
<< "AdjacencyEnergy "
<< "LaplacianTrace "
<< "LaplacianEnergy "
<< "LaplacianSpectralGap "
<< "NormalizedLaplacianTrace "
<< "NormalizedLaplacianEnergy "
<< "NormalizedLaplacianNumberOf2s "
<< "NormalizedLaplacianNumberOf1s "
<< "NormalizedLaplacianNumberOf0s "
<< "NormalizedLaplacianLowerSlope "
<< "NormalizedLaplacianUpperSlope "
<< "SmallWorldness"
<< std::endl;
std::stringstream localHeaderStream;
std::stringstream regionalHeaderStream;
std::stringstream globalDataStream;
std::stringstream localDataStream;
std::stringstream regionalDataStream;
std::string globalOutName = outName + "_global.txt";
std::string localOutName = outName + "_local.txt";
std::string regionalOutName = outName + "_regional.txt";
bool firstRun( true );
// iterate over all three possible methods
for(unsigned int method( 0 ); method < 3; method++)
{
// 0 - Random removal threshold
// 1 - Largest density below threshold
// 2 - Threshold based
// iterate over possible targets
for( unsigned int step( 0 ); step < granularity; step++ )
{
double targetValue( 0.0 );
bool newStep( true );
switch ( method )
{
case mitk::ConnectomicsNetworkThresholder::RandomRemovalOfWeakest :
case mitk::ConnectomicsNetworkThresholder::LargestLowerThanDensity :
targetValue = startDensity * (1 - static_cast<double>( step ) / ( granularity + 0.5 ) );
break;
case mitk::ConnectomicsNetworkThresholder::ThresholdBased :
targetValue = static_cast<double>( thresholdStepSize * step );
break;
default:
MITK_ERROR << "Invalid thresholding method called, aborting.";
return EXIT_FAILURE;
break;
}
mitk::ConnectomicsNetworkThresholder::Pointer thresholder = mitk::ConnectomicsNetworkThresholder::New();
thresholder->SetNetwork( network );
thresholder->SetTargetThreshold( targetValue );
thresholder->SetTargetDensity( targetValue );
thresholder->SetThresholdingScheme( static_cast<mitk::ConnectomicsNetworkThresholder::ThresholdingSchemes>(method) );
mitk::ConnectomicsNetwork::Pointer thresholdedNetwork = thresholder->GetThresholdedNetwork();
mitk::ConnectomicsStatisticsCalculator::Pointer statisticsCalculator = mitk::ConnectomicsStatisticsCalculator::New();
statisticsCalculator->SetNetwork( thresholdedNetwork );
statisticsCalculator->Update();
// global statistics
if( !noGlobalStatistics )
{
globalDataStream << statisticsCalculator->GetNumberOfVertices() << " "
<< statisticsCalculator->GetNumberOfEdges() << " "
<< statisticsCalculator->GetAverageDegree() << " "
<< statisticsCalculator->GetConnectionDensity() << " "
<< statisticsCalculator->GetNumberOfConnectedComponents() << " "
<< statisticsCalculator->GetAverageComponentSize() << " "
<< statisticsCalculator->GetLargestComponentSize() << " "
<< statisticsCalculator->GetRatioOfNodesInLargestComponent() << " "
<< statisticsCalculator->GetHopPlotExponent() << " "
<< statisticsCalculator->GetEffectiveHopDiameter() << " "
<< statisticsCalculator->GetAverageClusteringCoefficientsC() << " "
<< statisticsCalculator->GetAverageClusteringCoefficientsD() << " "
<< statisticsCalculator->GetAverageClusteringCoefficientsE() << " "
<< statisticsCalculator->GetAverageVertexBetweennessCentrality() << " "
<< statisticsCalculator->GetAverageEdgeBetweennessCentrality() << " "
<< statisticsCalculator->GetNumberOfIsolatedPoints() << " "
<< statisticsCalculator->GetRatioOfIsolatedPoints() << " "
<< statisticsCalculator->GetNumberOfEndPoints() << " "
<< statisticsCalculator->GetRatioOfEndPoints() << " "
<< statisticsCalculator->GetDiameter() << " "
<< statisticsCalculator->GetDiameter90() << " "
<< statisticsCalculator->GetRadius() << " "
<< statisticsCalculator->GetRadius90() << " "
<< statisticsCalculator->GetAverageEccentricity() << " "
<< statisticsCalculator->GetAverageEccentricity90() << " "
<< statisticsCalculator->GetAveragePathLength() << " "
<< statisticsCalculator->GetNumberOfCentralPoints() << " "
<< statisticsCalculator->GetRatioOfCentralPoints() << " "
<< statisticsCalculator->GetSpectralRadius() << " "
<< statisticsCalculator->GetSecondLargestEigenValue() << " "
<< statisticsCalculator->GetAdjacencyTrace() << " "
<< statisticsCalculator->GetAdjacencyEnergy() << " "
<< statisticsCalculator->GetLaplacianTrace() << " "
<< statisticsCalculator->GetLaplacianEnergy() << " "
<< statisticsCalculator->GetLaplacianSpectralGap() << " "
<< statisticsCalculator->GetNormalizedLaplacianTrace() << " "
<< statisticsCalculator->GetNormalizedLaplacianEnergy() << " "
<< statisticsCalculator->GetNormalizedLaplacianNumberOf2s() << " "
<< statisticsCalculator->GetNormalizedLaplacianNumberOf1s() << " "
<< statisticsCalculator->GetNormalizedLaplacianNumberOf0s() << " "
<< statisticsCalculator->GetNormalizedLaplacianLowerSlope() << " "
<< statisticsCalculator->GetNormalizedLaplacianUpperSlope() << " "
<< statisticsCalculator->GetSmallWorldness()
<< std::endl;
} // end global statistics
//create connectivity matrix png
if( createConnectivityMatriximage )
{
std::string connectivity_png_postfix = "_connectivity";
if( binaryConnectivity )
{
connectivity_png_postfix += "_binary";
}
else if( rescaleConnectivity )
{
connectivity_png_postfix += "_rescaled";
}
connectivity_png_postfix += ".png";
/* File format
* A png file depicting the binary connectivity matrix
*/
itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::Pointer filter = itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::New();
filter->SetInputNetwork( network );
filter->SetBinaryConnectivity( binaryConnectivity );
filter->SetRescaleConnectivity( rescaleConnectivity );
filter->Update();
typedef itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::OutputImageType connectivityMatrixImageType;
itk::ImageFileWriter< connectivityMatrixImageType >::Pointer connectivityWriter = itk::ImageFileWriter< connectivityMatrixImageType >::New();
connectivityWriter->SetInput( filter->GetOutput() );
connectivityWriter->SetFileName( outName + connectivity_png_postfix);
connectivityWriter->Update();
std::cout << "Connectivity matrix image written.";
} // end create connectivity matrix png
/*
* We can either calculate local indices for specific nodes, or specific regions
*/
// Create LabelToIndex translation
std::map< std::string, int > labelToIdMap;
std::vector< mitk::ConnectomicsNetwork::NetworkNode > nodeVector = thresholdedNetwork->GetVectorOfAllNodes();
for(int loop(0); loop < nodeVector.size(); loop++)
{
labelToIdMap.insert( std::pair< std::string, int>(nodeVector.at(loop).label, nodeVector.at(loop).id) );
}
std::vector< int > degreeVector = thresholdedNetwork->GetDegreeOfNodes();
std::vector< double > ccVector = thresholdedNetwork->GetLocalClusteringCoefficients( );
std::vector< double > bcVector = thresholdedNetwork->GetNodeBetweennessVector( );
// calculate local indices
{
// only add to header for the first step of the first method
if( firstRun )
{
localHeaderStream << "Th_method " << "Th_target " << "density";
}
double density = statisticsCalculator->GetConnectionDensity();
localDataStream << "\n"
<< method << " "
<< targetValue << " "
<< density;
for(unsigned int loop(0); loop < localLabels.size(); loop++ )
{
if( network->CheckForLabel(localLabels.at( loop )) )
{
if( firstRun )
{
localHeaderStream << " "
<< localLabels.at( loop ) << "_Degree "
<< localLabels.at( loop ) << "_CC "
<< localLabels.at( loop ) << "_BC";
}
localDataStream << " " << degreeVector.at( labelToIdMap.find( localLabels.at( loop ) )->second )
<< " " << ccVector.at( labelToIdMap.find( localLabels.at( loop ) )->second )
<< " " << bcVector.at( labelToIdMap.find( localLabels.at( loop ) )->second );
}
else
{
MITK_ERROR << "Illegal label. Label: \"" << localLabels.at( loop ) << "\" not found.";
}
}
}
// calculate regional indices
{
// only add to header for the first step of the first method
if( firstRun )
{
regionalHeaderStream << "Th_method " << "Th_target " << "density";
}
double density = statisticsCalculator->GetConnectionDensity();
regionalDataStream << "\n"
<< method << " "
<< targetValue << " "
<< density;
for( parsedRegionsIterator = parsedRegions.begin(); parsedRegionsIterator != parsedRegions.end(); parsedRegionsIterator++ )
{
std::vector<std::string> regionLabelsVector = parsedRegionsIterator->second;
std::string regionName = parsedRegionsIterator->first;
double sumDegree( 0 );
double sumCC( 0 );
double sumBC( 0 );
double count( 0 );
for( int loop(0); loop < regionLabelsVector.size(); loop++ )
{
if( thresholdedNetwork->CheckForLabel(regionLabelsVector.at( loop )) )
{
sumDegree = sumDegree + degreeVector.at( labelToIdMap.find( regionLabelsVector.at( loop ) )->second );
sumCC = sumCC + ccVector.at( labelToIdMap.find( regionLabelsVector.at( loop ) )->second );
sumBC = sumBC + bcVector.at( labelToIdMap.find( regionLabelsVector.at( loop ) )->second );
count = count + 1;
}
else
{
MITK_ERROR << "Illegal label. Label: \"" << regionLabelsVector.at( loop ) << "\" not found.";
}
}
// only add to header for the first step of the first method
if( firstRun )
{
regionalHeaderStream << " " << regionName << "_LocalAverageDegree "
<< regionName << "_LocalAverageCC "
<< regionName << "_LocalAverageBC "
<< regionName << "_NumberOfNodes";
}
regionalDataStream << " " << sumDegree / count
<< " " << sumCC / count
<< " " << sumBC / count
<< " " << count;
}
}
firstRun = false;
}
}// end calculate local averages
if( !noGlobalStatistics )
{
std::cout << "Writing to " << globalOutName;
std::ofstream glocalOutFile( globalOutName.c_str(), ios::out );
if( ! glocalOutFile.is_open() )
{
std::string errorMessage = "Could not open " + globalOutName + " for writing.";
MITK_ERROR << errorMessage;
return EXIT_FAILURE;
}
glocalOutFile << globalHeaderStream.str() << globalDataStream.str();
glocalOutFile.close();
}
if( localLabels.size() > 0 )
{
std::cout << "Writing to " << localOutName;
std::ofstream localOutFile( localOutName.c_str(), ios::out );
if( ! localOutFile.is_open() )
{
std::string errorMessage = "Could not open " + localOutName + " for writing.";
MITK_ERROR << errorMessage;
return EXIT_FAILURE;
}
localOutFile << localHeaderStream.str() << localDataStream.str();
localOutFile.close();
}
if( parsedRegions.size() > 0 )
{
std::cout << "Writing to " << regionalOutName;
std::ofstream regionalOutFile( regionalOutName.c_str(), ios::out );
if( ! regionalOutFile.is_open() )
{
std::string errorMessage = "Could not open " + regionalOutName + " for writing.";
MITK_ERROR << errorMessage;
return EXIT_FAILURE;
}
regionalOutFile << regionalHeaderStream.str() << regionalDataStream.str();
regionalOutFile.close();
}
return EXIT_SUCCESS;
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
std::cout << "DONE";
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp
index bebbbee478..7f8d7c47c0 100755
--- a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp
@@ -1,374 +1,374 @@
/*===================================================================
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 <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkFiniteDiffOdfMaximaExtractionFilter.h>
#include <mitkBaseDataIOFactory.h>
-#include <mitkDiffusionImage.h>
+#include <mitkImage.h>
#include <mitkQBallImage.h>
#include <mitkImageCast.h>
#include <mitkImageToItk.h>
#include <mitkTensorImage.h>
#include <mitkCoreObjectFactory.h>
#include "mitkCommandLineParser.h"
#include <itkShCoefficientImageImporter.h>
#include <itkFlipImageFilter.h>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include <mitkIOUtil.h>
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() )
{
std::cout << "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[])
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false);
parser.addArgument("outroot", "o", mitkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false);
parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask", "mask image");
parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true);
parser.addArgument("numpeaks", "p", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true);
parser.addArgument("peakthres", "r", mitkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true);
parser.addArgument("abspeakthres", "a", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true);
parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true);
parser.addArgument("noFlip", "f", mitkCommandLineParser::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"]);
std::cout << "image: " << imageName;
std::cout << "outroot: " << outRoot;
if (!maskImageName.empty())
std::cout << "mask: " << maskImageName;
else
std::cout << "no mask image selected";
std::cout << "numpeaks: " << numPeaks;
std::cout << "peakthres: " << peakThres;
std::cout << "abspeakthres: " << absPeakThres;
std::cout << "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;
std::cout << "Using FSL SH-basis";
}
else if ( boost::algorithm::equals(convention, "MRtrix") )
{
toolkitConvention = 2;
std::cout << "Using MRtrix SH-basis";
}
else
std::cout << "Using MITK SH-basis";
}
else
std::cout << "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)
{
std::cout << "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
{
std::cout << "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);
}
std::cout << "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(...)
{
std::cout << "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;
}
std::cout << "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");
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");
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::IOUtil::Save(directions.GetPointer(),outfilename.c_str());
}
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false);
parser.addArgument("shOrder", "sh", mitkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order");
parser.addArgument("outroot", "o", mitkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false);
parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask", "mask image");
parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true);
parser.addArgument("numpeaks", "p", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true);
parser.addArgument("peakthres", "r", mitkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true);
parser.addArgument("abspeakthres", "a", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true);
parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true);
parser.addArgument("noFlip", "f", mitkCommandLineParser::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;
}
diff --git a/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp b/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp
index 65cb7e0e33..13ea2b526b 100644
--- a/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp
@@ -1,239 +1,263 @@
/*===================================================================
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 "mitkBaseDataIOFactory.h"
#include <mitkCoreObjectFactory.h>
-#include "mitkDiffusionImage.h"
+#include "mitkImage.h"
#include "itkAnalyticalDiffusionQballReconstructionImageFilter.h"
#include <boost/lexical_cast.hpp>
#include "mitkCommandLineParser.h"
#include <mitkIOUtil.h>
#include <itksys/SystemTools.hxx>
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkITKImageImport.h>
+#include <mitkImageCast.h>
+#include <mitkProperties.h>
using namespace mitk;
/**
* Perform Q-ball reconstruction using a spherical harmonics basis
*/
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false);
parser.addArgument("outFile", "o", mitkCommandLineParser::OutputFile, "Output file", "output file", us::Any(), false);
parser.addArgument("shOrder", "sh", mitkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order", 4, true);
parser.addArgument("b0Threshold", "t", mitkCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0, true);
parser.addArgument("lambda", "r", mitkCommandLineParser::Float, "Lambda", "ragularization factor lambda", 0.006, true);
parser.addArgument("csa", "csa", mitkCommandLineParser::Bool, "Constant solid angle consideration", "use constant solid angle consideration");
parser.addArgument("outputCoeffs", "shc", mitkCommandLineParser::Bool, "Output coefficients", "output file containing the SH coefficients");
parser.addArgument("mrtrix", "mb", mitkCommandLineParser::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
{
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);
+ Image::Pointer dwi = dynamic_cast<Image*>(infile.at(0).GetPointer());
+ mitk::DiffusionPropertyHelper propertyHelper(dwi);
+ propertyHelper.AverageRedundantGradients(0.001);
+ propertyHelper.InitializeImage();
mitk::QBallImage::Pointer image = mitk::QBallImage::New();
mitk::Image::Pointer coeffsImage = mitk::Image::New();
std::cout << "SH order: " << shOrder;
std::cout << "lambda: " << lambda;
std::cout << "B0 threshold: " << threshold;
switch ( shOrder )
{
case 4:
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
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;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
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;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
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;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
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;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
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:
{
std::cout << "Supplied SH order not supported. Using default order of 4.";
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi));
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::IOUtil::SaveBaseData(image, outfilename);
if (outCoeffs)
mitk::IOUtil::SaveImage(coeffsImage, coeffout);
}
catch ( itk::ExceptionObject &err)
{
std::cout << "Exception: " << err;
}
catch ( std::exception err)
{
std::cout << "Exception: " << err.what();
}
catch ( ... )
{
std::cout << "Exception!";
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp b/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp
index 5f1b3a55db..ecbe39b53b 100644
--- a/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp
@@ -1,191 +1,190 @@
/*===================================================================
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 <mitkIOUtil.h>
#include "mitkImage.h"
#include <mitkImageCast.h>
#include "mitkITKImageImport.h"
#include <iostream>
#include <fstream>
#include <mitkTensorImage.h>
-#include <mitkDiffusionImage.h>
+#include <mitkImage.h>
+#include <mitkDiffusionPropertyHelper.h>
#include "itkTensorDerivedMeasurementsFilter.h"
#include "itkDiffusionTensor3DReconstructionImageFilter.h"
#include "mitkCommandLineParser.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 main(int argc, char* argv[])
{
std::cout << "TensorDerivedMapsExtraction";
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("help", "h", mitkCommandLineParser::String, "Help", "Show this help text");
parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input file", "input dwi file", us::Any(),false);
parser.addArgument("out", "o", mitkCommandLineParser::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";
std::cout << "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::Image::Pointer diffusionImage = mitk::IOUtil::LoadImage(inputFile);
+
+ if (diffusionImage.IsNull() || !mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(diffusionImage)) // 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++)
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainerCopy = mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New();
+ for( mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::ConstIterator it = mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImage)->Begin(); it != mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImage)->End(); it++)
{
gradientContainerCopy->push_back(it.Value());
}
- tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() );
- tensorReconstructionFilter->SetBValue(vols->GetReferenceBValue());
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(diffusionImage, itkVectorImagePointer);
+
+ tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, itkVectorImagePointer );
+ tensorReconstructionFilter->SetBValue( mitk::DiffusionPropertyHelper::GetReferenceBValue( diffusionImage ) );
tensorReconstructionFilter->SetThreshold(50);
tensorReconstructionFilter->Update();
typedef itk::Image<itk::DiffusionTensor3D<TTensorPixelType>, 3> TensorImageType;
TensorImageType::Pointer tensorImage = tensorReconstructionFilter->GetOutput();
- tensorImage->SetDirection( vols->GetVectorImage()->GetDirection() );
+ tensorImage->SetDirection( itkVectorImagePointer->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;
}
diff --git a/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp b/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp
index 2b7968a1fd..a1b066cc54 100644
--- a/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp
+++ b/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp
@@ -1,98 +1,103 @@
/*===================================================================
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 "mitkBaseDataIOFactory.h"
-#include "mitkDiffusionImage.h"
+#include "mitkImage.h"
+#include <mitkImageCast.h>
#include "mitkBaseData.h"
+#include <mitkDiffusionPropertyHelper.h>
#include <itkDiffusionTensor3DReconstructionImageFilter.h>
#include <itkDiffusionTensor3D.h>
#include <itkImageFileWriter.h>
#include <itkNrrdImageIO.h>
#include "mitkCommandLineParser.h"
#include <itksys/SystemTools.hxx>
using namespace mitk;
/**
* Convert files from one ending to the other
*/
int main(int argc, char* argv[])
{
std::cout << "TensorReconstruction";
mitkCommandLineParser parser;
parser.setArgumentPrefix("--", "-");
parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false);
parser.addArgument("outFile", "o", mitkCommandLineParser::OutputFile, "Output file", "output file", us::Any(), false);
parser.addArgument("b0Threshold", "t", mitkCommandLineParser::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
{
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());
+ Image::Pointer dwi = dynamic_cast<Image*>(infile.at(0).GetPointer());
+
+ mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, float > TensorReconstructionImageFilterType;
TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer );
+ filter->SetBValue( mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi ));
filter->SetThreshold(threshold);
filter->Update();
// Save tensor image
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)
{
std::cout << "Exception: " << err;
}
catch ( std::exception err)
{
std::cout << "Exception: " << err.what();
}
catch ( ... )
{
std::cout << "Exception!";
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/mitkFileFormatConverter.cpp b/Modules/DiffusionImaging/MiniApps/mitkFileFormatConverter.cpp
index 1cac8e8832..f758f91601 100755
--- a/Modules/DiffusionImaging/MiniApps/mitkFileFormatConverter.cpp
+++ b/Modules/DiffusionImaging/MiniApps/mitkFileFormatConverter.cpp
@@ -1,84 +1,80 @@
/*===================================================================
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 <mitkImageCast.h>
-#include <mitkDiffusionImage.h>
+#include <mitkImage.h>
#include <mitkBaseDataIOFactory.h>
#include <mitkIOUtil.h>
#include <mitkFiberBundleX.h>
#include "mitkCommandLineParser.h"
using namespace mitk;
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setTitle("Format Converter");
parser.setCategory("Fiber Tracking and Processing Methods");
parser.setDescription("");
parser.setContributor("MBI");
parser.setArgumentPrefix("--", "-");
parser.addArgument("in", "i", mitkCommandLineParser::InputFile, "Input:", "input file", us::Any(), false);
parser.addArgument("out", "o", mitkCommandLineParser::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
{
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::IOUtil::Save(dynamic_cast<DiffusionImage<short>*>(baseData.GetPointer()), outName.c_str());
- }
- else if ( dynamic_cast<Image*>(baseData.GetPointer()) )
+ if ( dynamic_cast<Image*>(baseData.GetPointer()) )
{
mitk::IOUtil::Save(dynamic_cast<Image*>(baseData.GetPointer()), outName.c_str());
}
else if ( dynamic_cast<FiberBundleX*>(baseData.GetPointer()) )
{
mitk::IOUtil::Save(dynamic_cast<FiberBundleX*>(baseData.GetPointer()) ,outName.c_str());
}
else
std::cout << "File type currently not supported!";
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/MiniApps/mitkRegistration.cpp b/Modules/DiffusionImaging/MiniApps/mitkRegistration.cpp
index d262915ba7..30077621d1 100644
--- a/Modules/DiffusionImaging/MiniApps/mitkRegistration.cpp
+++ b/Modules/DiffusionImaging/MiniApps/mitkRegistration.cpp
@@ -1,473 +1,457 @@
/*===================================================================
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.
===================================================================*/
// CTK
#include "mitkCommandLineParser.h"
#include <mitkIOUtil.h>
#include <mitkRegistrationWrapper.h>
#include <mitkImage.h>
#include <mitkImageCast.h>
#include <mitkITKImageImport.h>
-#include <mitkDiffusionImage.h>
#include <mitkImageTimeSelector.h>
// ITK
#include <itksys/SystemTools.hxx>
#include <itkDirectory.h>
#include "itkLinearInterpolateImageFunction.h"
#include "itkWindowedSincInterpolateImageFunction.h"
#include "itkIdentityTransform.h"
#include "itkResampleImageFilter.h"
typedef std::vector<std::string> FileListType;
typedef itk::Image<double, 3> InputImageType;
static mitk::Image::Pointer ExtractFirstTS(mitk::Image* image, std::string fileType)
{
if (fileType == ".dwi")
return image;
mitk::ImageTimeSelector::Pointer selector = mitk::ImageTimeSelector::New();
selector->SetInput(image);
selector->SetTimeNr(0);
selector->UpdateLargestPossibleRegion();
mitk::Image::Pointer img =selector->GetOutput()->Clone();
return img;
}
static std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems)
{
std::stringstream ss(s);
std::string item;
while (std::getline(ss, item, delim))
{
elems.push_back(item);
}
return elems;
}
static std::vector<std::string> split(const std::string &s, char delim)
{
std::vector < std::string > elems;
return split(s, delim, elems);
}
/// Create list of all files in provided folder ending with same postfix
static FileListType CreateFileList(std::string folder , std::string postfix)
{
itk::Directory::Pointer dir = itk::Directory::New();
FileListType fileList;
if( dir->Load(folder.c_str() ) )
{
int n = dir->GetNumberOfFiles();
for(int r=0;r<n;r++)
{
std::string filename = dir->GetFile( r );
if (filename == "." || filename == "..")
continue;
filename = folder + filename;
if (!itksys::SystemTools::FileExists( filename.c_str()))
continue;
if (filename.substr(filename.length() -postfix.length() ) == postfix)
fileList.push_back(filename);
}
}
return fileList;
}
static std::string GetSavePath(std::string outputFolder, std::string fileName)
{
std::string fileType = itksys::SystemTools::GetFilenameExtension(fileName);
std::string fileStem = itksys::SystemTools::GetFilenameWithoutExtension(fileName);
std::string savePathAndFileName = outputFolder +fileStem + fileType;
return savePathAndFileName;
}
static mitk::Image::Pointer ResampleBySpacing(mitk::Image *input, float *spacing, bool useLinInt = false)
{
InputImageType::Pointer itkImage = InputImageType::New();
CastToItkImage(input,itkImage);
/**
* 1) Resampling
*
*/
// Identity transform.
// We don't want any transform on our image except rescaling which is not
// specified by a transform but by the input/output spacing as we will see
// later.
// So no transform will be specified.
typedef itk::IdentityTransform<double, 3> T_Transform;
// The resampler type itself.
typedef itk::ResampleImageFilter<InputImageType, InputImageType> T_ResampleFilter;
// Prepare the resampler.
// Instantiate the transform and specify it should be the id transform.
T_Transform::Pointer _pTransform = T_Transform::New();
_pTransform->SetIdentity();
// Instantiate the resampler. Wire in the transform and the interpolator.
T_ResampleFilter::Pointer _pResizeFilter = T_ResampleFilter::New();
_pResizeFilter->SetTransform(_pTransform);
// Set the output origin.
_pResizeFilter->SetOutputOrigin(itkImage->GetOrigin());
// Compute the size of the output.
// The size (# of pixels) in the output is recomputed using
// the ratio of the input and output sizes.
InputImageType::SpacingType inputSpacing = itkImage->GetSpacing();
InputImageType::SpacingType outputSpacing;
const InputImageType::RegionType& inputSize = itkImage->GetLargestPossibleRegion();
InputImageType::SizeType outputSize;
typedef InputImageType::SizeType::SizeValueType SizeValueType;
// Set the output spacing.
outputSpacing[0] = spacing[0];
outputSpacing[1] = spacing[1];
outputSpacing[2] = spacing[2];
outputSize[0] = static_cast<SizeValueType>(inputSize.GetSize()[0] * inputSpacing[0] / outputSpacing[0] + .5);
outputSize[1] = static_cast<SizeValueType>(inputSize.GetSize()[1] * inputSpacing[1] / outputSpacing[1] + .5);
outputSize[2] = static_cast<SizeValueType>(inputSize.GetSize()[2] * inputSpacing[2] / outputSpacing[2] + .5);
_pResizeFilter->SetOutputSpacing(outputSpacing);
_pResizeFilter->SetSize(outputSize);
typedef itk::LinearInterpolateImageFunction< InputImageType > LinearInterpolatorType;
LinearInterpolatorType::Pointer lin_interpolator = LinearInterpolatorType::New();
typedef itk::Function::WelchWindowFunction<4> WelchWindowFunction;
typedef itk::WindowedSincInterpolateImageFunction< InputImageType, 4,WelchWindowFunction> WindowedSincInterpolatorType;
WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New();
if (useLinInt)
_pResizeFilter->SetInterpolator(lin_interpolator);
else
_pResizeFilter->SetInterpolator(sinc_interpolator);
// Specify the input.
_pResizeFilter->SetInput(itkImage);
_pResizeFilter->Update();
mitk::Image::Pointer image = mitk::Image::New();
image->InitializeByItk(_pResizeFilter->GetOutput());
mitk::GrabItkImageMemory( _pResizeFilter->GetOutput(), image);
return image;
}
/// Build a derived file name from moving images e.g. xxx_T2.nrrd becomes xxx_GTV.nrrd
static FileListType CreateDerivedFileList(std::string baseFN, std::string baseSuffix, std::vector<std::string> derivedPatterns)
{
FileListType files;
for (unsigned int i=0; i < derivedPatterns.size(); i++)
{
std::string derResourceSuffix = derivedPatterns.at(i);
std::string derivedResourceFilename = baseFN.substr(0,baseFN.length() -baseSuffix.length()) + derResourceSuffix;
MITK_INFO <<" Looking for file: " << derivedResourceFilename;
if (!itksys::SystemTools::FileExists(derivedResourceFilename.c_str()))
{
MITK_INFO << "CreateDerivedFileList: File does not exit. Skipping entry.";
continue;
}
files.push_back(derivedResourceFilename);
}
return files;
}
/// Save images according to file type
static void SaveImage(std::string fileName, mitk::Image* image, std::string fileType )
{
MITK_INFO << "----Save to " << fileName;
- if (fileType == "dwi") // IOUtil does not handle dwi files properly Bug 15772
- {
- try
- {
- mitk::IOUtil::Save(dynamic_cast<mitk::DiffusionImage<short>*>(image), fileName.c_str());
- }
- catch( const itk::ExceptionObject& e)
- {
- MITK_ERROR << "Caught exception: " << e.what();
- mitkThrow() << "Failed with exception from subprocess!";
- }
- }
- else
- {
- mitk::IOUtil::SaveImage(image, fileName);
- }
+ mitk::IOUtil::Save(image, fileName);
}
/// Copy derived resources from first time step. Append _reg tag, but leave data untouched.
static void CopyResources(FileListType fileList, std::string outputPath)
{
for (unsigned int j=0; j < fileList.size(); j++)
{
std::string derivedResourceFilename = fileList.at(j);
std::string fileType = itksys::SystemTools::GetFilenameExtension(derivedResourceFilename);
std::string fileStem = itksys::SystemTools::GetFilenameWithoutExtension(derivedResourceFilename);
std::string savePathAndFileName = outputPath +fileStem + "." + fileType;
MITK_INFO << "Copy resource " << savePathAndFileName;
mitk::Image::Pointer resImage = ExtractFirstTS(mitk::IOUtil::LoadImage(derivedResourceFilename), fileType);
mitk::IOUtil::SaveImage(resImage, savePathAndFileName);
}
}
int main( int argc, char* argv[] )
{
mitkCommandLineParser parser;
parser.setArgumentPrefix("--","-");
parser.setTitle("Folder Registraton");
parser.setCategory("Preprocessing Tools");
parser.setDescription("http://docs.mitk.org/nightly-qt4/DiffusionMiniApps.html");
parser.setContributor("MBI");
// Add command line argument names
parser.addArgument("help", "h",mitkCommandLineParser::Bool, "Help", "Show this help text");
//parser.addArgument("usemask", "u", QVariant::Bool, "Use segmentations (derived resources) to exclude areas from registration metrics");
parser.addArgument("input", "i", mitkCommandLineParser::InputDirectory, "Input:", "Input folder",us::Any(),false);
parser.addArgument("output", "o", mitkCommandLineParser::OutputDirectory, "Output:", "Output folder (ending with /)",us::Any(),false);
parser.addArgument("fixed", "f", mitkCommandLineParser::String, "Fixed images:", "Suffix for fixed image (if none is supplied first file matching moving pattern is chosen)",us::Any(),true);
parser.addArgument("moving", "m", mitkCommandLineParser::String, "Moving images:", "Suffix for moving images",us::Any(),false);
parser.addArgument("derived", "d", mitkCommandLineParser::String, "Derived resources:", "Derived resources suffixes (replaces suffix for moving images); comma separated",us::Any(),true);
parser.addArgument("silent", "s", mitkCommandLineParser::Bool, "Silent:" "No xml progress output.");
parser.addArgument("resample", "r", mitkCommandLineParser::String, "Resample (x,y,z)mm:", "Resample provide x,y,z spacing in mm (e.g. -r 1,1,3), is not applied to tensor data",us::Any());
parser.addArgument("binary", "b", mitkCommandLineParser::Bool, "Binary:", "Speficies that derived resource are binary (interpolation using nearest neighbor)",us::Any());
parser.addArgument("correct-origin", "c", mitkCommandLineParser::Bool, "Origin correction:", "Correct for large origin displacement. Switch when you reveive: Joint PDF summed to zero ",us::Any());
parser.addArgument("sinc-int", "s", mitkCommandLineParser::Bool, "Windowed-sinc interpolation:", "Use windowed-sinc interpolation (3) instead of linear interpolation ",us::Any());
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
// Handle special arguments
bool silent = false;
bool isBinary = false;
bool alignOrigin = false;
bool useLinearInterpol = true;
{
if (parsedArgs.size() == 0)
{
return EXIT_FAILURE;
}
if (parsedArgs.count("sinc-int"))
useLinearInterpol = false;
if (parsedArgs.count("silent"))
silent = true;
if (parsedArgs.count("binary"))
isBinary = true;
if (parsedArgs.count("correct-origin"))
alignOrigin = true;
// Show a help message
if ( parsedArgs.count("help") || parsedArgs.count("h"))
{
std::cout << parser.helpText();
return EXIT_SUCCESS;
}
}
std::string refPattern = "";
bool useFirstMoving = false;
std::string movingImgPattern = us::any_cast<string>(parsedArgs["moving"]);
if (parsedArgs.count("fixed"))
{
refPattern = us::any_cast<string>(parsedArgs["fixed"]);
}
else
{
useFirstMoving = true;
refPattern = movingImgPattern;
}
std::string outputPath = us::any_cast<string>(parsedArgs["output"]);
std::string inputPath = us::any_cast<string>(parsedArgs["input"]);
//QString resampleReference = parsedArgs["resample"].toString();
//bool maskTumor = parsedArgs["usemask"].toBool();
// if derived sources pattern is provided, populate QStringList with possible filename postfixes
std::vector<std::string> derPatterns;
if (parsedArgs.count("derived") || parsedArgs.count("d") )
{
std::string arg = us::any_cast<string>(parsedArgs["derived"]);
derPatterns = split(arg ,',');
}
std::vector<std::string> spacings;
float spacing[3];
bool doResampling = false;
if (parsedArgs.count("resample") || parsedArgs.count("d") )
{
std::string arg = us::any_cast<string>(parsedArgs["resample"]);
spacings = split(arg ,',');
spacing[0] = atoi(spacings.at(0).c_str());
spacing[1] = atoi(spacings.at(1).c_str());
spacing[2] = atoi(spacings.at(2).c_str());
doResampling = true;
}
MITK_INFO << "Input Folder : " << inputPath;
MITK_INFO << "Looking for reference image ...";
FileListType referenceFileList = CreateFileList(inputPath,refPattern);
if ((!useFirstMoving && referenceFileList.size() != 1) || (useFirstMoving && referenceFileList.size() == 0))
{
MITK_ERROR << "None or more than one possible reference images (" << refPattern <<") found. Exiting." << referenceFileList.size();
MITK_INFO << "Choose a fixed arguement that is unique in the given folder!";
return EXIT_FAILURE;
}
std::string referenceFileName = referenceFileList.at(0);
MITK_INFO << "Loading Reference (fixed) image: " << referenceFileName;
std::string fileType = itksys::SystemTools::GetFilenameExtension(referenceFileName);
mitk::Image::Pointer refImage = ExtractFirstTS(mitk::IOUtil::LoadImage(referenceFileName), fileType);
mitk::Image::Pointer resampleReference = NULL;
if (doResampling)
{
refImage = ResampleBySpacing(refImage,spacing);
resampleReference = refImage;
}
if (refImage.IsNull())
MITK_ERROR << "Loaded fixed image is NULL";
// Copy reference image to destination
std::string savePathAndFileName = GetSavePath(outputPath, referenceFileName);
mitk::IOUtil::SaveImage(refImage, savePathAndFileName);
// Copy all derived resources also to output folder, adding _reg suffix
referenceFileList = CreateDerivedFileList(referenceFileName, movingImgPattern,derPatterns);
CopyResources(referenceFileList, outputPath);
std::string derivedResourceFilename;
mitk::Image::Pointer referenceMask = NULL; // union of all segmentations
if (!silent)
{
// XML Output to report progress
std::cout << "<filter-start>";
std::cout << "<filter-name>Batched Registration</filter-name>";
std::cout << "<filter-comment>Starting registration ... </filter-comment>";
std::cout << "</filter-start>";
}
// Now iterate over all files and register them to the reference image,
// also register derived resources based on file patterns
// ------------------------------------------------------------------------------
// Create File list
FileListType movingImagesList = CreateFileList(inputPath, movingImgPattern);
// TODO Reactivate Resampling Feature
// mitk::Image::Pointer resampleImage = NULL;
// if (QFileInfo(resampleReference).isFile())
// {
// resampleImage = mitk::IOUtil::LoadImage(resampleReference.toStdString());
// }
for (unsigned int i =0; i < movingImagesList.size(); i++)
{
std::string fileMorphName = movingImagesList.at(i);
if (fileMorphName == referenceFileName)
{
// do not process reference image again
continue;
}
MITK_INFO << "Processing image " << fileMorphName;
// 1 Register morphological file to reference image
if (!itksys::SystemTools::FileExists(fileMorphName.c_str()))
{
MITK_WARN << "File does not exit. Skipping entry.";
continue;
}
// Origin of images is cancelled
// TODO make this optional!!
double transf[6];
double offset[3];
{
std::string fileType = itksys::SystemTools::GetFilenameExtension(fileMorphName);
mitk::Image::Pointer movingImage = ExtractFirstTS(mitk::IOUtil::LoadImage(fileMorphName), fileType);
if (movingImage.IsNull())
MITK_ERROR << "Loaded moving image is NULL";
// Store transformation, apply it to morph file
MITK_INFO << "----------Registering moving image to reference----------";
mitk::RegistrationWrapper::GetTransformation(refImage, movingImage, transf, offset, alignOrigin, referenceMask);
mitk::RegistrationWrapper::ApplyTransformationToImage(movingImage, transf,offset, resampleReference); // , resampleImage
savePathAndFileName = GetSavePath(outputPath, fileMorphName);
if (fileType == ".dwi")
fileType = "dwi";
SaveImage(savePathAndFileName,movingImage,fileType );
}
if (!silent)
{
std::cout << "<filter-progress-text progress=\"" <<
(float)i / (float)movingImagesList.size()
<< "\" >.</filter-progress-text>";
}
// Now parse all derived resource and apply the above calculated transformation to them
// ------------------------------------------------------------------------------------
FileListType fList = CreateDerivedFileList(fileMorphName, movingImgPattern,derPatterns);
if (fList.size() > 0)
MITK_INFO << "----------DERIVED RESOURCES ---------";
for (unsigned int j=0; j < fList.size(); j++)
{
derivedResourceFilename = fList.at(j);
MITK_INFO << "----Processing derived resorce " << derivedResourceFilename << " ...";
std::string fileType = itksys::SystemTools::GetFilenameExtension(derivedResourceFilename);
mitk::Image::Pointer derivedMovingResource = ExtractFirstTS(mitk::IOUtil::LoadImage(derivedResourceFilename), fileType);
// Apply transformation to derived resource, treat derived resource as binary
mitk::RegistrationWrapper::ApplyTransformationToImage(derivedMovingResource, transf,offset, resampleReference,isBinary);
savePathAndFileName = GetSavePath(outputPath, derivedResourceFilename);
SaveImage(savePathAndFileName,derivedMovingResource,fileType );
}
}
if (!silent)
std::cout << "<filter-end/>";
return EXIT_SUCCESS;
}