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