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