diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberBundleXReaderWriterTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberBundleXReaderWriterTest.cpp index ca8c7b2a91..6e048ee7fe 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberBundleXReaderWriterTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberBundleXReaderWriterTest.cpp @@ -1,69 +1,71 @@ /*=================================================================== 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 "mitkTestingMacros.h" #include #include #include #include #include #include "mitkTestFixture.h" class mitkFiberBundleXReaderWriterTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkFiberBundleXReaderWriterTestSuite); MITK_TEST(Equal_SaveLoad_ReturnsTrue); CPPUNIT_TEST_SUITE_END(); private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ mitk::FiberBundleX::Pointer fib1; mitk::FiberBundleX::Pointer fib2; public: void setUp() { fib1 = NULL; fib2 = NULL; std::string filename = GetTestDataFilePath("DiffusionImaging/fiberBundleX.fib"); - mitk::BaseData::Pointer baseData = mitk::IOUtil::LoadBaseData(filename); + std::vector fibInfile = mitk::IOUtil::Load( filename); + mitk::BaseData::Pointer baseData = fibInfile.at(0); + fib1 = dynamic_cast(baseData.GetPointer()); } void tearDown() { fib1 = NULL; fib2 = NULL; } void Equal_SaveLoad_ReturnsTrue() { mitk::IOUtil::Save(fib1.GetPointer(), std::string(MITK_TEST_OUTPUT_DIR)+"/writerTest.fib"); - mitk::BaseData::Pointer baseData = mitk::IOUtil::LoadBaseData(std::string(MITK_TEST_OUTPUT_DIR)+"/writerTest.fib"); - fib2 = dynamic_cast(baseData.GetPointer()); + std::vector baseData = mitk::IOUtil::Load(std::string(MITK_TEST_OUTPUT_DIR)+"/writerTest.fib"); + fib2 = dynamic_cast(baseData[0].GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", fib1->Equals(fib2)); //MITK_ASSERT_EQUAL(fib1, fib2, "A saved and re-loaded file should be equal"); } }; MITK_TEST_SUITE_REGISTRATION(mitkFiberBundleXReaderWriter) diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkGibbsTrackingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkGibbsTrackingTest.cpp index cc8d8c0c10..cbfe9626ed 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkGibbsTrackingTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkGibbsTrackingTest.cpp @@ -1,92 +1,92 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include +#include using namespace mitk; /**Documentation * Test for gibbs tracking filter */ int mitkGibbsTrackingTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkGibbsTrackingTest"); MITK_TEST_CONDITION_REQUIRED(argc>4,"check for input data") QBallImage::Pointer mitkQballImage; Image::Pointer mitkMaskImage; mitk::FiberBundleX::Pointer fib1; try{ MITK_INFO << "Q-Ball image: " << argv[1]; MITK_INFO << "Mask image: " << argv[2]; MITK_INFO << "Parameter file: " << argv[3]; MITK_INFO << "Reference bundle: " << argv[4]; - const std::string s1="", s2=""; - std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[1], s1, s2, false ); + std::vector infile = mitk::IOUtil::Load( argv[1]); mitkQballImage = dynamic_cast(infile.at(0).GetPointer()); MITK_TEST_CONDITION_REQUIRED(mitkQballImage.IsNotNull(),"check qball image") - infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[2], s1, s2, false ); + infile = mitk::IOUtil::Load( argv[2] ); mitkMaskImage = dynamic_cast(infile.at(0).GetPointer()); MITK_TEST_CONDITION_REQUIRED(mitkMaskImage.IsNotNull(),"check mask image") - infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[4], s1, s2, false ); + infile = mitk::IOUtil::Load( argv[4]); fib1 = dynamic_cast(infile.at(0).GetPointer()); MITK_TEST_CONDITION_REQUIRED(fib1.IsNotNull(),"check fiber bundle") typedef itk::Vector OdfVectorType; typedef itk::Image OdfVectorImgType; typedef itk::Image MaskImgType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; OdfVectorImgType::Pointer itk_qbi = OdfVectorImgType::New(); mitk::CastToItkImage(mitkQballImage, itk_qbi); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); gibbsTracker->SetQBallImage(itk_qbi.GetPointer()); gibbsTracker->SetMaskImage(itk_mask); gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetRandomSeed(1); gibbsTracker->SetLoadParameterFile(argv[3]); gibbsTracker->Update(); mitk::FiberBundleX::Pointer fib2 = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); MITK_TEST_CONDITION_REQUIRED(fib1->Equals(fib2), "check if gibbs tracking has changed"); gibbsTracker->SetRandomSeed(0); gibbsTracker->Update(); fib2 = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); MITK_TEST_CONDITION_REQUIRED(!fib1->Equals(fib2), "check if gibbs tracking has changed after wrong seed"); } catch(...) { return EXIT_FAILURE; } // always end with this! MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkPeakExtractionTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkPeakExtractionTest.cpp index 6110dc612f..6e1f252e20 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkPeakExtractionTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkPeakExtractionTest.cpp @@ -1,108 +1,106 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include -#include #include #include #include #include #include #include #include #include #include #include #include using namespace std; int mitkPeakExtractionTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkStreamlineTrackingTest"); MITK_TEST_CONDITION_REQUIRED(argc>3,"check for input data") string shCoeffFileName = argv[1]; string maskFileName = argv[2]; string referenceFileName = argv[3]; MITK_INFO << "SH-coefficient file: " << shCoeffFileName; MITK_INFO << "Mask file: " << maskFileName; MITK_INFO << "Reference fiber file: " << referenceFileName; try { mitk::CoreObjectFactory::GetInstance(); mitk::Image::Pointer image = mitk::IOUtil::LoadImage(shCoeffFileName); mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(maskFileName); typedef itk::Image ItkUcharImgType; typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, 4, 20242 > MaximaExtractionFilterType; MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); MITK_INFO << "Casting mask image ..."; ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); mitk::CastToItkImage(mitkMaskImage, itkMask); filter->SetMaskImage(itkMask); MITK_INFO << "Casting SH image ..."; typedef mitk::ImageToItk< MaximaExtractionFilterType::CoefficientImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); filter->SetInput(caster->GetOutput()); filter->SetMaxNumPeaks(2); filter->SetPeakThreshold(0.4); filter->SetAbsolutePeakThreshold(0.01); filter->SetAngularThreshold(25); filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); filter->SetNumberOfThreads(1); MITK_INFO << "Starting extraction ..."; filter->Update(); mitk::FiberBundleX::Pointer fib1 = filter->GetOutputFiberBundle(); MITK_INFO << "Loading reference ..."; - const std::string s1="", s2=""; - std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( referenceFileName, s1, s2, false ); + std::vector infile = mitk::IOUtil::Load( referenceFileName ); mitk::FiberBundleX::Pointer fib2 = dynamic_cast(infile.at(0).GetPointer()); // TODO: reduce epsilon. strange issues with differing values between windows and linux. MITK_TEST_CONDITION_REQUIRED(fib1->Equals(fib2), "Check if tractograms are equal."); } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp b/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp index fb7e4fa586..75dce2cdd5 100644 --- a/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp +++ b/Modules/DiffusionImaging/MiniApps/DiffusionIndices.cpp @@ -1,143 +1,143 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include /** * Calculate indices derived from Qball or tensor images */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Diffusion Indices"); parser.setCategory("Diffusion Related Measures"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", us::Any(), false); parser.addArgument("index", "idx", mitkCommandLineParser::String, "Index:", "index (fa, gfa, ra, ad, rd, ca, l2, l3, md)", us::Any(), false); parser.addArgument("outFile", "o", mitkCommandLineParser::OutputFile, "Output:", "output file", us::Any(), false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string inFileName = us::any_cast(parsedArgs["input"]); string index = us::any_cast(parsedArgs["index"]); string outFileName = us::any_cast(parsedArgs["outFile"]); string ext = itksys::SystemTools::GetFilenameLastExtension(outFileName); if (ext.empty()) outFileName += ".nrrd"; try { // load input image - mitk::BaseData::Pointer infile = mitk::IOUtil::LoadBaseData( inFileName ); + std::vector infile = mitk::IOUtil::Load( inFileName ); if( boost::algorithm::ends_with(inFileName, ".qbi") && index=="gfa" ) { typedef itk::Vector OdfVectorType; typedef itk::Image ItkQballImageType; - mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(infile.GetPointer()); + mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(infile[0].GetPointer()); ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New(); mitk::CastToItkImage(mitkQballImage, itk_qbi); typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(itk_qbi); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(gfaFilter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if( boost::algorithm::ends_with(inFileName, ".dti") ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; - mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.GetPointer()); + mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile[0].GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itk_dti.GetPointer() ); if(index=="fa") measurementsCalculator->SetMeasure(MeasurementsType::FA); else if(index=="ra") measurementsCalculator->SetMeasure(MeasurementsType::RA); else if(index=="ad") measurementsCalculator->SetMeasure(MeasurementsType::AD); else if(index=="rd") measurementsCalculator->SetMeasure(MeasurementsType::RD); else if(index=="ca") measurementsCalculator->SetMeasure(MeasurementsType::CA); else if(index=="l2") measurementsCalculator->SetMeasure(MeasurementsType::L2); else if(index=="l3") measurementsCalculator->SetMeasure(MeasurementsType::L3); else if(index=="md") measurementsCalculator->SetMeasure(MeasurementsType::MD); else { MITK_WARN << "No valid diffusion index for input image (tensor image) defined"; return EXIT_FAILURE; } measurementsCalculator->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(measurementsCalculator->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else std::cout << "Diffusion index " << index << " not supported for supplied file type."; } 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/FiberProcessing.cpp b/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp index 5bbb301414..dd839bb3b2 100644 --- a/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp +++ b/Modules/DiffusionImaging/MiniApps/FiberProcessing.cpp @@ -1,206 +1,204 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include mitk::FiberBundleX::Pointer LoadFib(std::string filename) { - const std::string s1="", s2=""; - std::vector fibInfile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false ); + std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; - mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } int main(int argc, char* argv[]) { std::cout << "FiberProcessing"; mitkCommandLineParser parser; parser.setTitle("Fiber Processing"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib)", us::Any(), false); parser.addArgument("outFile", "o", mitkCommandLineParser::OutputFile, "Output:", "output fiber bundle (.fib)", us::Any(), false); parser.addArgument("smooth", "s", mitkCommandLineParser::Float, "Spline resampling:", "Resample fiber using splines with the given point distance (in mm)"); parser.addArgument("compress", "c", mitkCommandLineParser::Float, "Compress:", "Compress fiber using the given error threshold (in mm)"); parser.addArgument("minLength", "l", mitkCommandLineParser::Float, "Minimum length:", "Minimum fiber length (in mm)"); parser.addArgument("maxLength", "m", mitkCommandLineParser::Float, "Maximum length:", "Maximum fiber length (in mm)"); parser.addArgument("minCurv", "a", mitkCommandLineParser::Float, "Minimum curvature radius:", "Minimum curvature radius (in mm)"); parser.addArgument("mirror", "p", mitkCommandLineParser::Int, "Invert coordinates:", "Invert fiber coordinates XYZ (e.g. 010 to invert y-coordinate of each fiber point)"); parser.addArgument("rotate-x", "rx", mitkCommandLineParser::Float, "Rotate x-axis:", "Rotate around x-axis (if copy is given the copy is rotated, in deg)"); parser.addArgument("rotate-y", "ry", mitkCommandLineParser::Float, "Rotate y-axis:", "Rotate around y-axis (if copy is given the copy is rotated, in deg)"); parser.addArgument("rotate-z", "rz", mitkCommandLineParser::Float, "Rotate z-axis:", "Rotate around z-axis (if copy is given the copy is rotated, in deg)"); parser.addArgument("scale-x", "sx", mitkCommandLineParser::Float, "Scale x-axis:", "Scale in direction of x-axis (if copy is given the copy is scaled)"); parser.addArgument("scale-y", "sy", mitkCommandLineParser::Float, "Scale y-axis:", "Scale in direction of y-axis (if copy is given the copy is scaled)"); parser.addArgument("scale-z", "sz", mitkCommandLineParser::Float, "Scale z-axis", "Scale in direction of z-axis (if copy is given the copy is scaled)"); parser.addArgument("translate-x", "tx", mitkCommandLineParser::Float, "Translate x-axis:", "Translate in direction of x-axis (if copy is given the copy is translated, in mm)"); parser.addArgument("translate-y", "ty", mitkCommandLineParser::Float, "Translate y-axis:", "Translate in direction of y-axis (if copy is given the copy is translated, in mm)"); parser.addArgument("translate-z", "tz", mitkCommandLineParser::Float, "Translate z-axis:", "Translate in direction of z-axis (if copy is given the copy is translated, in mm)"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; float smoothDist = -1; if (parsedArgs.count("smooth")) smoothDist = us::any_cast(parsedArgs["smooth"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); float minFiberLength = -1; if (parsedArgs.count("minLength")) minFiberLength = us::any_cast(parsedArgs["minLength"]); float maxFiberLength = -1; if (parsedArgs.count("maxLength")) maxFiberLength = us::any_cast(parsedArgs["maxLength"]); float curvThres = -1; if (parsedArgs.count("minCurv")) curvThres = us::any_cast(parsedArgs["minCurv"]); int axis = 0; if (parsedArgs.count("mirror")) axis = us::any_cast(parsedArgs["mirror"]); float rotateX = 0; if (parsedArgs.count("rotate-x")) rotateX = us::any_cast(parsedArgs["rotate-x"]); float rotateY = 0; if (parsedArgs.count("rotate-y")) rotateY = us::any_cast(parsedArgs["rotate-y"]); float rotateZ = 0; if (parsedArgs.count("rotate-z")) rotateZ = us::any_cast(parsedArgs["rotate-z"]); float scaleX = 0; if (parsedArgs.count("scale-x")) scaleX = us::any_cast(parsedArgs["scale-x"]); float scaleY = 0; if (parsedArgs.count("scale-y")) scaleY = us::any_cast(parsedArgs["scale-y"]); float scaleZ = 0; if (parsedArgs.count("scale-z")) scaleZ = us::any_cast(parsedArgs["scale-z"]); float translateX = 0; if (parsedArgs.count("translate-x")) translateX = us::any_cast(parsedArgs["translate-x"]); float translateY = 0; if (parsedArgs.count("translate-y")) translateY = us::any_cast(parsedArgs["translate-y"]); float translateZ = 0; if (parsedArgs.count("translate-z")) translateZ = us::any_cast(parsedArgs["translate-z"]); string inFileName = us::any_cast(parsedArgs["input"]); string outFileName = us::any_cast(parsedArgs["outFile"]); try { mitk::FiberBundleX::Pointer fib = LoadFib(inFileName); if (minFiberLength>0) fib->RemoveShortFibers(minFiberLength); if (maxFiberLength>0) fib->RemoveLongFibers(maxFiberLength); if (curvThres>0) fib->ApplyCurvatureThreshold(curvThres, false); if (smoothDist>0) fib->ResampleSpline(smoothDist); if (compress>0) fib->Compress(compress); if (axis/100==1) fib->MirrorFibers(0); if ((axis%100)/10==1) fib->MirrorFibers(1); if (axis%10==1) fib->MirrorFibers(2); if (rotateX > 0 || rotateY > 0 || rotateZ > 0){ std::cout << "Rotate " << rotateX << " " << rotateY << " " << rotateZ; fib->RotateAroundAxis(rotateX, rotateY, rotateZ); } if (translateX > 0 || translateY > 0 || translateZ > 0){ fib->TranslateFibers(translateX, translateY, translateZ); } if (scaleX > 0 || scaleY > 0 || scaleZ > 0) fib->ScaleFibers(scaleX, scaleY, scaleZ); mitk::IOUtil::SaveBaseData(fib.GetPointer(), outFileName ); } 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/FileFormatConverter.cpp b/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp index 5725c33e21..e9ae47d575 100644 --- a/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp +++ b/Modules/DiffusionImaging/MiniApps/FileFormatConverter.cpp @@ -1,79 +1,77 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include "mitkCommandLineParser.h" 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 parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string inName = us::any_cast(parsedArgs["in"]); string outName = us::any_cast(parsedArgs["out"]); try { - const std::string s1="", s2=""; - std::vector infile = BaseDataIO::LoadBaseDataFromFile( inName, s1, s2, false ); - mitk::BaseData::Pointer baseData = infile.at(0); + std::vector baseData = mitk::IOUtil::Load(inName); - if ( dynamic_cast(baseData.GetPointer()) ) + if ( baseData.size()>0 && dynamic_cast(baseData[0].GetPointer()) ) { - mitk::IOUtil::Save(dynamic_cast(baseData.GetPointer()), outName.c_str()); + mitk::IOUtil::Save(dynamic_cast(baseData[0].GetPointer()), outName.c_str()); } - else if ( dynamic_cast(baseData.GetPointer()) ) + else if ( baseData.size()>0 && dynamic_cast(baseData[0].GetPointer()) ) { - mitk::IOUtil::Save(dynamic_cast(baseData.GetPointer()) ,outName.c_str()); + mitk::IOUtil::Save(dynamic_cast(baseData[0].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/GibbsTracking.cpp b/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp index 07ab216380..eb159cae12 100755 --- a/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp +++ b/Modules/DiffusionImaging/MiniApps/GibbsTracking.cpp @@ -1,241 +1,238 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include template typename itk::ShCoefficientImageImporter< float, shOrder >::QballImageType::Pointer TemplatedConvertShCoeffs(mitk::Image* mitkImg, int toolkit, bool noFlip = false) { typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImg); caster->Update(); itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput(); typename FilterType::Pointer filter = FilterType::New(); if (noFlip) { filter->SetInputImage(itkImage); } else { std::cout << "Flipping image"; itk::FixedArray flipAxes; flipAxes[0] = true; flipAxes[1] = true; flipAxes[2] = false; flipAxes[3] = false; itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New(); flipper->SetInput(itkImage); flipper->SetFlipAxes(flipAxes); flipper->Update(); itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput(); itk::Matrix< double,4,4 > m = itkImage->GetDirection(); m[0][0] *= -1; m[1][1] *= -1; flipped->SetDirection(m); itk::Point< float, 4 > o = itkImage->GetOrigin(); o[0] -= (flipped->GetLargestPossibleRegion().GetSize(0)-1); o[1] -= (flipped->GetLargestPossibleRegion().GetSize(1)-1); flipped->SetOrigin(o); filter->SetInputImage(flipped); } switch (toolkit) { case 0: filter->SetToolkit(FilterType::FSL); break; case 1: filter->SetToolkit(FilterType::MRTRIX); break; default: filter->SetToolkit(FilterType::FSL); } filter->GenerateData(); return filter->GetQballImage(); } int main(int argc, char* argv[]) { std::cout << "GibbsTracking"; mitkCommandLineParser parser; parser.setTitle("Gibbs Tracking"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", us::Any(), false); parser.addArgument("parameters", "p", mitkCommandLineParser::InputFile, "Parameters:", "parameter file (.gtp)", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "binary mask image"); parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "SH coefficient:", "sh coefficient convention (FSL, MRtrix)", string("FSL"), true); parser.addArgument("outFile", "o", mitkCommandLineParser::OutputFile, "Output:", "output fiber bundle (.fib)", us::Any(), false); parser.addArgument("noFlip", "f", mitkCommandLineParser::Bool, "No flip:", "do not flip input image to match MITK coordinate convention"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string inFileName = us::any_cast(parsedArgs["input"]); string paramFileName = us::any_cast(parsedArgs["parameters"]); string outFileName = us::any_cast(parsedArgs["outFile"]); bool noFlip = false; if (parsedArgs.count("noFlip")) noFlip = us::any_cast(parsedArgs["noFlip"]); try { // instantiate gibbs tracker typedef itk::Vector OdfVectorType; typedef itk::Image ItkQballImageType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); // load input image - const std::string s1="", s2=""; - std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false ); - - mitk::Image::Pointer mitkImage = dynamic_cast(infile.at(0).GetPointer()); + mitk::Image::Pointer mitkImage = mitk::IOUtil::LoadImage(inFileName); // try to cast to qball image if( boost::algorithm::ends_with(inFileName, ".qbi") ) { std::cout << "Loading qball image ..."; - mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(infile.at(0).GetPointer()); + mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(mitkImage.GetPointer()); ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New(); mitk::CastToItkImage(mitkQballImage, itk_qbi); gibbsTracker->SetQBallImage(itk_qbi.GetPointer()); } else if( boost::algorithm::ends_with(inFileName, ".dti") ) { std::cout << "Loading tensor image ..."; typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; - mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); + mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(mitkImage.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); gibbsTracker->SetTensorImage(itk_dti); } else if ( boost::algorithm::ends_with(inFileName, ".nii") ) { std::cout << "Loading sh-coefficient image ..."; int nrCoeffs = mitkImage->GetLargestPossibleRegion().GetSize()[3]; int c=3, d=2-2*nrCoeffs; double D = c*c-4*d; int shOrder; if (D>0) { shOrder = (-c+sqrt(D))/2.0; if (shOrder<0) shOrder = (-c-sqrt(D))/2.0; } else if (D==0) shOrder = -c/2.0; std::cout << "using SH-order " << shOrder; int toolkitConvention = 0; if (parsedArgs.count("shConvention")) { string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); if ( boost::algorithm::equals(convention, "MRtrix") ) { toolkitConvention = 1; std::cout << "Using MRtrix style sh-coefficient convention"; } else std::cout << "Using FSL style sh-coefficient convention"; } else std::cout << "Using FSL style sh-coefficient convention"; switch (shOrder) { case 4: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<4>(mitkImage, toolkitConvention, noFlip)); break; case 6: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<6>(mitkImage, toolkitConvention, noFlip)); break; case 8: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<8>(mitkImage, toolkitConvention, noFlip)); break; case 10: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<10>(mitkImage, toolkitConvention, noFlip)); break; case 12: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<12>(mitkImage, toolkitConvention, noFlip)); break; default: std::cout << "SH-order " << shOrder << " not supported"; } } else return EXIT_FAILURE; // global tracking if (parsedArgs.count("mask")) { typedef itk::Image MaskImgType; mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(us::any_cast(parsedArgs["mask"])); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); gibbsTracker->SetMaskImage(itk_mask); } gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetLoadParameterFile( paramFileName ); // gibbsTracker->SetLutPath( "" ); gibbsTracker->Update(); mitk::FiberBundleX::Pointer mitkFiberBundle = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); mitkFiberBundle->SetReferenceGeometry(mitkImage->GetGeometry()); mitk::IOUtil::SaveBaseData(mitkFiberBundle.GetPointer(), outFileName ); } 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/ImageResampler.cpp b/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp index a05ce2419b..5c49782177 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageResampler.cpp @@ -1,314 +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 #include #include #include #include #include #include #include // ITK #include #include #include "itkLinearInterpolateImageFunction.h" #include "itkWindowedSincInterpolateImageFunction.h" #include "itkIdentityTransform.h" #include "itkResampleImageFilter.h" #include "itkResampleDwiImageFilter.h" typedef itk::Image InputImageType; 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 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 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 &split(const std::string &s, char delim, std::vector &elems) { std::stringstream ss(s); std::string item; while (std::getline(ss, item, delim)) { elems.push_back(item); } return elems; } static std::vector 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 T_Transform; // The resampler type itself. typedef itk::ResampleImageFilter 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(inputSize.GetSize()[0] * inputSpacing[0] / outputSpacing[0] + .5); outputSize[1] = static_cast(inputSize.GetSize()[1] * inputSpacing[1] / outputSpacing[1] + .5); outputSize[2] = static_cast(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; mitk::IOUtil::Save(image, fileName); } mitk::Image::Pointer ResampleDWIbySpacing(mitk::Image::Pointer input, float* spacing, bool useLinInt = true) { itk::Vector spacingVector; spacingVector[0] = spacing[0]; spacingVector[1] = spacing[1]; spacingVector[2] = spacing[2]; typedef itk::ResampleDwiImageFilter ResampleFilterType; mitk::DiffusionPropertyHelper::ImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::ImageType::New(); mitk::CastToItkImage(input, itkVectorImagePointer); ResampleFilterType::Pointer resampler = ResampleFilterType::New(); resampler->SetInput( itkVectorImagePointer ); resampler->SetInterpolation(ResampleFilterType::Interpolate_Linear); resampler->SetNewSpacing(spacingVector); resampler->Update(); 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 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(parsedArgs["output"]); std::string inputFile = us::any_cast(parsedArgs["input"]); std::vector spacings; float spacing[3]; if (parsedArgs.count("spacing")) { std::string arg = us::any_cast(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(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); - mitk::Image::Pointer inputDWI = dynamic_cast(mitk::IOUtil::LoadBaseData(inputFile).GetPointer()); - if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(inputDWI)) + mitk::Image::Pointer inputDWI = mitk::IOUtil::LoadImage(inputFile); + if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(inputDWI.GetPointer())) { 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 fa556c06b5..0b99f56cae 100644 --- a/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp +++ b/Modules/DiffusionImaging/MiniApps/MultishellMethods.cpp @@ -1,218 +1,216 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include 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 parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string inName = us::any_cast(parsedArgs["in"]); string outName = us::any_cast(parsedArgs["out"]); bool applyADC = us::any_cast(parsedArgs["adc"]); bool applyAKC = us::any_cast(parsedArgs["akc"]); bool applyBiExp = us::any_cast(parsedArgs["biexp"]); string targetType = us::any_cast(parsedArgs["targetbvalue"]); try { std::cout << "Loading " << inName; - const std::string s1="", s2=""; - std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( inName, s1, s2, false ); - mitk::BaseData::Pointer baseData = infile.at(0); - if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(baseData.GetPointer()) ) ) + mitk::Image::Pointer dwi = mitk::IOUtil::LoadImage(inName); + + if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dwi ) ) { - mitk::Image::Pointer dwi = dynamic_cast(baseData.GetPointer()); typedef itk::RadialMultishellToSingleshellImageFilter FilterType; typedef itk::DwiGradientLengthCorrectionFilter CorrectionFilterType; CorrectionFilterType::Pointer roundfilter = CorrectionFilterType::New(); roundfilter->SetRoundingValue( 1000 ); roundfilter->SetReferenceBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi )); roundfilter->SetReferenceGradientDirectionContainer(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi)); roundfilter->Update(); 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::DiffusionPropertyHelper::BValueMapType &originalShellMap = mitk::DiffusionPropertyHelper::GetBValueMap(dwi); mitk::DiffusionPropertyHelper::ImageType::Pointer vectorImage = mitk::DiffusionPropertyHelper::ImageType::New(); mitk::CastToItkImage(dwi, vectorImage); const mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainer = mitk::DiffusionPropertyHelper::GetGradientContainer(dwi); const unsigned int &bValue = mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi ); // filter call vnl_vector bValueList(originalShellMap.size()-1); double targetBValue = bValueList.mean(); 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::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::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::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/NetworkCreation.cpp b/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp index d824593643..751aff65ca 100644 --- a/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp +++ b/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp @@ -1,132 +1,132 @@ /*=================================================================== 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 // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include "mitkConnectomicsNetworkCreator.h" #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("fiberImage", "f", mitkCommandLineParser::InputFile, "Input image", "input fiber image (.fib)", us::Any(), false); parser.addArgument("parcellation", "p", mitkCommandLineParser::InputFile, "Parcellation image", "parcellation image", us::Any(), false); parser.addArgument("outputNetwork", "o", mitkCommandLineParser::String, "Output network", "where to save the output (.cnf)", us::Any(), false); parser.addArgument("radius", "r", mitkCommandLineParser::Int, "Radius", "Search radius in mm", 15, true); parser.addArgument("noCenterOfMass", "com", mitkCommandLineParser::Bool, "No center of mass", "Do not use center of mass for node positions"); parser.setCategory("Connectomics"); parser.setTitle("Network Creation"); parser.setDescription(""); parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; //default values int searchRadius( 15 ); bool noCenterOfMass( false ); // parse command line arguments std::string fiberFilename = us::any_cast(parsedArgs["fiberImage"]); std::string parcellationFilename = us::any_cast(parsedArgs["parcellation"]); std::string outputFilename = us::any_cast(parsedArgs["outputNetwork"]); if (parsedArgs.count("radius")) searchRadius = us::any_cast(parsedArgs["radius"]); if (parsedArgs.count("noCenterOfMass")) noCenterOfMass = us::any_cast(parsedArgs["noCenterOfMass"]); try { const std::string s1="", s2=""; // load fiber image std::vector fiberInfile = - mitk::BaseDataIO::LoadBaseDataFromFile( fiberFilename, s1, s2, false ); + mitk::IOUtil::Load( fiberFilename); if( fiberInfile.empty() ) { std::string errorMessage = "Fiber Image at " + fiberFilename + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* fiberBaseData = fiberInfile.at(0); mitk::FiberBundleX* fiberBundle = dynamic_cast( fiberBaseData ); // load parcellation std::vector parcellationInFile = - mitk::BaseDataIO::LoadBaseDataFromFile( parcellationFilename, s1, s2, false ); + mitk::IOUtil::Load( parcellationFilename); if( parcellationInFile.empty() ) { std::string errorMessage = "Parcellation at " + parcellationFilename + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* parcellationBaseData = parcellationInFile.at(0); mitk::Image* parcellationImage = dynamic_cast( parcellationBaseData ); // do creation mitk::ConnectomicsNetworkCreator::Pointer connectomicsNetworkCreator = mitk::ConnectomicsNetworkCreator::New(); connectomicsNetworkCreator->SetSegmentation( parcellationImage ); connectomicsNetworkCreator->SetFiberBundle( fiberBundle ); if( !noCenterOfMass ) { connectomicsNetworkCreator->CalculateCenterOfMass(); } connectomicsNetworkCreator->SetEndPointSearchRadius( searchRadius ); connectomicsNetworkCreator->CreateNetworkFromFibersAndSegmentation(); mitk::ConnectomicsNetwork::Pointer network = connectomicsNetworkCreator->GetNetwork(); std::cout << "searching writer"; mitk::IOUtil::SaveBaseData(network.GetPointer(), outputFilename ); 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/NetworkStatistics.cpp b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp index d4e43e56f3..1abe56c0bf 100644 --- a/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp +++ b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp @@ -1,515 +1,514 @@ /*=================================================================== 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 #include #include #include #include #include // boost includes #include // ITK includes #include // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include #include #include +#include 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 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(parsedArgs["inputNetwork"]); std::string outName = us::any_cast(parsedArgs["outputFile"]); mitkCommandLineParser::StringContainerType localLabels; if(parsedArgs.count("localStatistics")) { localLabels = us::any_cast(parsedArgs["localStatistics"]); } mitkCommandLineParser::StringContainerType unparsedRegions; std::map< std::string, std::vector > parsedRegions; std::map< std::string, std::vector >::iterator parsedRegionsIterator; if(parsedArgs.count("regionList")) { unparsedRegions = us::any_cast(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 >( tempRegionVector.at(0), insertRegionVector) ); } else { MITK_ERROR << "Region already exists. Skipping second occurrence."; } } } if (parsedArgs.count("noGlobalStatistics")) noGlobalStatistics = us::any_cast(parsedArgs["noGlobalStatistics"]); if (parsedArgs.count("binaryConnectivity")) binaryConnectivity = us::any_cast(parsedArgs["binaryConnectivity"]); if (parsedArgs.count("rescaleConnectivity")) rescaleConnectivity = us::any_cast(parsedArgs["rescaleConnectivity"]); if (parsedArgs.count("createConnectivityMatriximage")) createConnectivityMatriximage = us::any_cast(parsedArgs["createConnectivityMatriximage"]); if (parsedArgs.count("granularity")) granularity = us::any_cast(parsedArgs["granularity"]); if (parsedArgs.count("startDensity")) startDensity = us::any_cast(parsedArgs["startDensity"]); if (parsedArgs.count("thresholdStepSize")) thresholdStepSize = us::any_cast(parsedArgs["thresholdStepSize"]); try { - const std::string s1="", s2=""; - // load network std::vector networkFile = - mitk::BaseDataIO::LoadBaseDataFromFile( networkName, s1, s2, false ); + mitk::IOUtil::Load( networkName); 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( 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( step ) / ( granularity + 0.5 ) ); break; case mitk::ConnectomicsNetworkThresholder::ThresholdBased : targetValue = static_cast( 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(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 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 7f8d7c47c0..fcec556ffc 100755 --- a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp +++ b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp @@ -1,374 +1,355 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include -#include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include -mitk::Image::Pointer LoadData(std::string filename) -{ - if( filename.empty() ) - return NULL; - - const std::string s1="", s2=""; - std::vector 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(baseData.GetPointer()); -} - - template 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 parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string imageName = us::any_cast(parsedArgs["image"]); string outRoot = us::any_cast(parsedArgs["outroot"]); // optional arguments string maskImageName(""); if (parsedArgs.count("mask")) maskImageName = us::any_cast(parsedArgs["mask"]); int normalization = 1; if (parsedArgs.count("normalization")) normalization = us::any_cast(parsedArgs["normalization"]); int numPeaks = 2; if (parsedArgs.count("numpeaks")) numPeaks = us::any_cast(parsedArgs["numpeaks"]); float peakThres = 0.4; if (parsedArgs.count("peakthres")) peakThres = us::any_cast(parsedArgs["peakthres"]); float absPeakThres = 0.06; if (parsedArgs.count("abspeakthres")) absPeakThres = us::any_cast(parsedArgs["abspeakthres"]); bool noFlip = false; if (parsedArgs.count("noFlip")) noFlip = us::any_cast(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); + mitk::Image::Pointer image = mitk::IOUtil::LoadImage(imageName); + mitk::Image::Pointer mask = mitk::IOUtil::LoadImage(maskImageName); typedef itk::Image 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(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 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; iSize(); 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(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 parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; int shOrder = -1; if (parsedArgs.count("shOrder")) shOrder = us::any_cast(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; }