diff --git a/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp b/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp index 0e36fe7952..394cb74058 100644 --- a/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp +++ b/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp @@ -1,139 +1,143 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" // std includes #include // CTK includes #include "ctkCommandLineParser.h" // MITK includes #include #include "mitkConnectomicsNetworkCreator.h" #include int NetworkCreation(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("fiberImage", "f", ctkCommandLineParser::String, "input fiber image (.fib)", us::Any(), false); - parser.addArgument("parcellation", "p", ctkCommandLineParser::String, "parcellation image", us::Any(), false); - parser.addArgument("outputNetwork", "o", ctkCommandLineParser::String, "where to save the ouput (.cnf)", us::Any(), false); + parser.addArgument("fiberImage", "f", ctkCommandLineParser::File, "Input image", "input fiber image (.fib)", us::Any(), false); + parser.addArgument("parcellation", "p", ctkCommandLineParser::File, "Parcellation image", "parcellation image", us::Any(), false); + parser.addArgument("outputNetwork", "o", ctkCommandLineParser::String, "Output network", "where to save the output (.cnf)", us::Any(), false); - parser.addArgument("radius", "r", ctkCommandLineParser::Int, "Search radius in mm", 15, true); - parser.addArgument("noCenterOfMass", "com", ctkCommandLineParser::Bool, "Do not use center of mass for node positions"); + parser.addArgument("radius", "r", ctkCommandLineParser::Int, "Radius", "Search radius in mm", 15, true); + parser.addArgument("noCenterOfMass", "com", ctkCommandLineParser::Bool, "No center of mass", "Do not use center of mass for node positions"); + parser.setCategory("Connectomics"); + parser.setTitle("Network Creation"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; //default values int searchRadius( 15 ); bool noCenterOfMass( false ); // parse command line arguments std::string fiberFilename = us::any_cast(parsedArgs["fiberImage"]); std::string parcellationFilename = us::any_cast(parsedArgs["parcellation"]); std::string outputFilename = us::any_cast(parsedArgs["outputNetwork"]); if (parsedArgs.count("radius")) searchRadius = us::any_cast(parsedArgs["radius"]); if (parsedArgs.count("noCenterOfMass")) noCenterOfMass = us::any_cast(parsedArgs["noCenterOfMass"]); try { const std::string s1="", s2=""; // load fiber image std::vector fiberInfile = mitk::BaseDataIO::LoadBaseDataFromFile( fiberFilename, s1, s2, false ); if( fiberInfile.empty() ) { std::string errorMessage = "Fiber Image at " + fiberFilename + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* fiberBaseData = fiberInfile.at(0); mitk::FiberBundleX* fiberBundle = dynamic_cast( fiberBaseData ); // load parcellation std::vector parcellationInFile = mitk::BaseDataIO::LoadBaseDataFromFile( parcellationFilename, s1, s2, false ); if( parcellationInFile.empty() ) { std::string errorMessage = "Parcellation at " + parcellationFilename + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* parcellationBaseData = parcellationInFile.at(0); mitk::Image* parcellationImage = dynamic_cast( parcellationBaseData ); // do creation mitk::ConnectomicsNetworkCreator::Pointer connectomicsNetworkCreator = mitk::ConnectomicsNetworkCreator::New(); connectomicsNetworkCreator->SetSegmentation( parcellationImage ); connectomicsNetworkCreator->SetFiberBundle( fiberBundle ); if( !noCenterOfMass ) { connectomicsNetworkCreator->CalculateCenterOfMass(); } connectomicsNetworkCreator->SetEndPointSearchRadius( searchRadius ); connectomicsNetworkCreator->CreateNetworkFromFibersAndSegmentation(); mitk::ConnectomicsNetwork::Pointer network = connectomicsNetworkCreator->GetNetwork(); MITK_INFO << "searching writer"; mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) { if ( (*it)->CanWriteBaseDataType(network.GetPointer()) ) { MITK_INFO << "writing"; (*it)->SetFileName( outputFilename.c_str() ); (*it)->DoWrite( network.GetPointer() ); } } return EXIT_SUCCESS; } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_INFO << "DONE"; return EXIT_SUCCESS; } RegisterDiffusionMiniApp(NetworkCreation); diff --git a/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp index 8b8163b9b7..e395e6e6ec 100644 --- a/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp +++ b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp @@ -1,513 +1,518 @@ /*=================================================================== 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 // std includes #include #include #include #include #include #include // boost includes #include // ITK includes #include // CTK includes #include "ctkCommandLineParser.h" // MITK includes #include #include #include #include int NetworkStatistics(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("inputNetwork", "i", ctkCommandLineParser::String, "input connectomics network (.cnf)", us::Any(), false); - parser.addArgument("outputFile", "o", ctkCommandLineParser::String, "name of output file", us::Any(), false); - - parser.addArgument("noGlobalStatistics", "g", ctkCommandLineParser::Bool, "Do not calculate global statistics"); - parser.addArgument("createConnectivityMatriximage", "I", ctkCommandLineParser::Bool, "Write connectivity matrix image"); - parser.addArgument("binaryConnectivity", "b", ctkCommandLineParser::Bool, "Whether to create a binary connectivity matrix"); - parser.addArgument("rescaleConnectivity", "r", ctkCommandLineParser::Bool, "Whether to rescale the connectivity matrix"); - parser.addArgument("localStatistics", "L", ctkCommandLineParser::StringList, "Provide a list of node labels for local statistics", us::Any()); - parser.addArgument("regionList", "R", ctkCommandLineParser::StringList, "A space separated list of regions. Each region has the format\n regionname;label1;label2;...;labelN", us::Any()); - parser.addArgument("granularity", "gr", ctkCommandLineParser::Int, "How finely to test the density range and how many thresholds to consider"); - parser.addArgument("startDensity", "d", ctkCommandLineParser::Bool, "Largest density for the range"); - parser.addArgument("thresholdStepSize", "t", ctkCommandLineParser::Int, "Distance of two adjacent thresholds"); + parser.addArgument("inputNetwork", "i", ctkCommandLineParser::File, "Input network", "input connectomics network (.cnf)", us::Any(), false); + parser.addArgument("outputFile", "o", ctkCommandLineParser::File, "Output file", "name of output file", us::Any(), false); + + parser.addArgument("noGlobalStatistics", "g", ctkCommandLineParser::Bool, "No global statistics", "Do not calculate global statistics"); + parser.addArgument("createConnectivityMatriximage", "I", ctkCommandLineParser::Bool, "Write connectivity matrix image", "Write connectivity matrix image"); + parser.addArgument("binaryConnectivity", "b", ctkCommandLineParser::Bool, "Binary connectivity", "Whether to create a binary connectivity matrix"); + parser.addArgument("rescaleConnectivity", "r", ctkCommandLineParser::Bool, "Rescale connectivity", "Whether to rescale the connectivity matrix"); + parser.addArgument("localStatistics", "L", ctkCommandLineParser::StringList, "Local statistics", "Provide a list of node labels for local statistics", us::Any()); + parser.addArgument("regionList", "R", ctkCommandLineParser::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", ctkCommandLineParser::Int, "Granularity", "How finely to test the density range and how many thresholds to consider"); + parser.addArgument("startDensity", "d", ctkCommandLineParser::Bool, "Start Density", "Largest density for the range"); + parser.addArgument("thresholdStepSize", "t", ctkCommandLineParser::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 ); unsigned int granularity( 1 ); double startDensity( 1.0 ); unsigned int thresholdStepSize( 3 ); // parse command line arguments std::string networkName = us::any_cast(parsedArgs["inputNetwork"]); std::string outName = us::any_cast(parsedArgs["outputFile"]); ctkCommandLineParser::StringContainerType localLabels; if(parsedArgs.count("localStatistics")) { localLabels = us::any_cast(parsedArgs["localStatistics"]); } ctkCommandLineParser::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(); MITK_INFO << "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 ) { MITK_INFO << "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 ) { MITK_INFO << "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 ) { MITK_INFO << "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) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_INFO << "DONE"; return EXIT_SUCCESS; } RegisterDiffusionMiniApp(NetworkStatistics); diff --git a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp index 4d4b92e740..a19d6a791c 100755 --- a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp +++ b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp @@ -1,372 +1,377 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include #include #include #include #include #include #include #include #include #include #include "ctkCommandLineParser.h" #include #include #include #include #include mitk::Image::Pointer LoadData(std::string filename) { if( filename.empty() ) return NULL; const std::string s1="", s2=""; std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false ); if( infile.empty() ) { MITK_INFO << "File " << filename << " could not be read!"; return NULL; } mitk::BaseData::Pointer baseData = infile.at(0); return dynamic_cast(baseData.GetPointer()); } template int StartPeakExtraction(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("image", "i", ctkCommandLineParser::String, "sh coefficient image", us::Any(), false); - parser.addArgument("outroot", "o", ctkCommandLineParser::String, "output root", us::Any(), false); - parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); - parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "0=no norm, 1=max norm, 2=single vec norm", 1, true); - parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "maximum number of extracted peaks", 2, true); - parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "peak threshold relative to largest peak", 0.4, true); - parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "absolute peak threshold weighted with local GFA value", 0.06, true); - parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); - parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "do not flip input image to match MITK coordinate convention"); + parser.addArgument("image", "i", ctkCommandLineParser::File, "Input image", "sh coefficient image", us::Any(), false); + parser.addArgument("outroot", "o", ctkCommandLineParser::Directory, "Output directory", "output root", us::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::File, "Mask", "mask image"); + parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); + parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); + parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true); + parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true); + parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); + parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention"); + + parser.setCategory("Preprocessing Tools"); + parser.setTitle("Peak Extraction"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string imageName = us::any_cast(parsedArgs["image"]); string outRoot = us::any_cast(parsedArgs["outroot"]); // optional arguments string maskImageName(""); if (parsedArgs.count("mask")) maskImageName = us::any_cast(parsedArgs["mask"]); int normalization = 1; if (parsedArgs.count("normalization")) normalization = us::any_cast(parsedArgs["normalization"]); int numPeaks = 2; if (parsedArgs.count("numpeaks")) numPeaks = us::any_cast(parsedArgs["numpeaks"]); float peakThres = 0.4; if (parsedArgs.count("peakthres")) peakThres = us::any_cast(parsedArgs["peakthres"]); float absPeakThres = 0.06; if (parsedArgs.count("abspeakthres")) absPeakThres = us::any_cast(parsedArgs["abspeakthres"]); bool noFlip = false; if (parsedArgs.count("noFlip")) noFlip = us::any_cast(parsedArgs["noFlip"]); MITK_INFO << "image: " << imageName; MITK_INFO << "outroot: " << outRoot; if (!maskImageName.empty()) MITK_INFO << "mask: " << maskImageName; else MITK_INFO << "no mask image selected"; MITK_INFO << "numpeaks: " << numPeaks; MITK_INFO << "peakthres: " << peakThres; MITK_INFO << "abspeakthres: " << absPeakThres; MITK_INFO << "shOrder: " << shOrder; try { mitk::Image::Pointer image = LoadData(imageName); mitk::Image::Pointer mask = LoadData(maskImageName); typedef itk::Image ItkUcharImgType; typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); int toolkitConvention = 0; if (parsedArgs.count("shConvention")) { string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); if ( boost::algorithm::equals(convention, "FSL") ) { toolkitConvention = 1; MITK_INFO << "Using FSL SH-basis"; } else if ( boost::algorithm::equals(convention, "MRtrix") ) { toolkitConvention = 2; MITK_INFO << "Using MRtrix SH-basis"; } else MITK_INFO << "Using MITK SH-basis"; } else MITK_INFO << "Using MITK SH-basis"; ItkUcharImgType::Pointer itkMaskImage = NULL; if (mask.IsNotNull()) { try{ itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(mask, itkMaskImage); filter->SetMaskImage(itkMaskImage); } catch(...) { } } if (toolkitConvention>0) { MITK_INFO << "Converting coefficient image to MITK format"; typedef itk::ShCoefficientImageImporter< float, shOrder > ConverterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput(); typename ConverterType::Pointer converter = ConverterType::New(); if (noFlip) { converter->SetInputImage(itkImage); } else { MITK_INFO << "Flipping image"; itk::FixedArray flipAxes; flipAxes[0] = true; flipAxes[1] = true; flipAxes[2] = false; flipAxes[3] = false; itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New(); flipper->SetInput(itkImage); flipper->SetFlipAxes(flipAxes); flipper->Update(); itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput(); itk::Matrix< double,4,4 > m = itkImage->GetDirection(); m[0][0] *= -1; m[1][1] *= -1; flipped->SetDirection(m); itk::Point< float, 4 > o = itkImage->GetOrigin(); o[0] -= (flipped->GetLargestPossibleRegion().GetSize(0)-1); o[1] -= (flipped->GetLargestPossibleRegion().GetSize(1)-1); flipped->SetOrigin(o); converter->SetInputImage(flipped); } MITK_INFO << "Starting conversion"; switch (toolkitConvention) { case 1: converter->SetToolkit(ConverterType::FSL); filter->SetToolkit(MaximaExtractionFilterType::FSL); break; case 2: converter->SetToolkit(ConverterType::MRTRIX); filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); break; default: converter->SetToolkit(ConverterType::FSL); filter->SetToolkit(MaximaExtractionFilterType::FSL); break; } converter->GenerateData(); filter->SetInput(converter->GetCoefficientImage()); } else { try{ typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); filter->SetInput(caster->GetOutput()); } catch(...) { MITK_INFO << "wrong image type"; return EXIT_FAILURE; } } filter->SetMaxNumPeaks(numPeaks); filter->SetPeakThreshold(peakThres); filter->SetAbsolutePeakThreshold(absPeakThres); filter->SetAngularThreshold(1); switch (normalization) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } MITK_INFO << "Starting extraction"; filter->Update(); // write direction images { typedef typename MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer; typename ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer(); for (unsigned int i=0; iSize(); i++) { typename MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i); if (itkMaskImage.IsNotNull()) { itkImg->SetDirection(itkMaskImage->GetDirection()); itkImg->SetOrigin(itkMaskImage->GetOrigin()); } string outfilename = outRoot; outfilename.append("_DIRECTION_"); outfilename.append(boost::lexical_cast(i)); outfilename.append(".nrrd"); MITK_INFO << "writing " << outfilename; typedef itk::ImageFileWriter< typename MaximaExtractionFilterType::ItkDirectionImage > WriterType; typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(itkImg); writer->Update(); } } // write num directions image { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); if (itkMaskImage.IsNotNull()) { numDirImage->SetDirection(itkMaskImage->GetDirection()); numDirImage->SetOrigin(itkMaskImage->GetOrigin()); } string outfilename = outRoot.c_str(); outfilename.append("_NUM_DIRECTIONS.nrrd"); MITK_INFO << "writing " << outfilename; typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(numDirImage); writer->Update(); } // write vector field { mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle(); string outfilename = outRoot.c_str(); outfilename.append("_VECTOR_FIELD.fib"); mitk::FiberBundleXWriter::Pointer fibWriter = mitk::FiberBundleXWriter::New(); fibWriter->SetFileName(outfilename.c_str()); fibWriter->DoWrite(directions.GetPointer()); } } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_INFO << "DONE"; return EXIT_SUCCESS; } int PeakExtraction(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", ctkCommandLineParser::String, "sh coefficient image", us::Any(), false); parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "spherical harmonics order"); parser.addArgument("outroot", "o", ctkCommandLineParser::String, "output root", us::Any(), false); parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "0=no norm, 1=max norm, 2=single vec norm", 1, true); parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "maximum number of extracted peaks", 2, true); parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "peak threshold relative to largest peak", 0.4, true); parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "absolute peak threshold weighted with local GFA value", 0.06, true); parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "do not flip input image to match MITK coordinate convention"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; int shOrder = -1; if (parsedArgs.count("shOrder")) shOrder = us::any_cast(parsedArgs["shOrder"]); switch (shOrder) { case 4: return StartPeakExtraction<4>(argc, argv); case 6: return StartPeakExtraction<6>(argc, argv); case 8: return StartPeakExtraction<8>(argc, argv); case 10: return StartPeakExtraction<10>(argc, argv); case 12: return StartPeakExtraction<12>(argc, argv); } return EXIT_FAILURE; } RegisterDiffusionMiniApp(PeakExtraction); diff --git a/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp b/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp index 5c000230a4..ef4f84b569 100755 --- a/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp +++ b/Modules/DiffusionImaging/MiniApps/PeaksAngularError.cpp @@ -1,205 +1,210 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include #include #include #include #include #include #include "ctkCommandLineParser.h" #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include int PeaksAngularError(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("test", "t", ctkCommandLineParser::StringList, "test direction images", us::Any(), false); - parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", us::Any(), false); - parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", us::Any(), false); - parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); - parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results"); - parser.addArgument("ignore", "i", ctkCommandLineParser::Bool, "don't increase error for missing or too many directions"); + parser.addArgument("test", "t", ctkCommandLineParser::StringList, "Test images", "test direction images", us::Any(), false); + parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "Reference images", "reference direction images", us::Any(), false); + parser.addArgument("out", "o", ctkCommandLineParser::Directory, "Output directory", "output root", us::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::File, "Mask", "mask image"); + parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "Verbose", "output optional and intermediate calculation results"); + parser.addArgument("ignore", "i", ctkCommandLineParser::Bool, "Ignore", "don't increase error for missing or too many directions"); + + parser.setCategory("Preprocessing Tools"); + parser.setTitle("Peaks Angular Error"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; ctkCommandLineParser::StringContainerType testImages = us::any_cast(parsedArgs["test"]); ctkCommandLineParser::StringContainerType referenceImages = us::any_cast(parsedArgs["reference"]); string maskImage(""); if (parsedArgs.count("mask")) maskImage = us::any_cast(parsedArgs["mask"]); string outRoot = us::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); bool ignore = false; if (parsedArgs.count("ignore")) ignore = us::any_cast(parsedArgs["ignore"]); try { typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; ItkDirectionImageContainerType::Pointer directionImageContainer = ItkDirectionImageContainerType::New(); for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(testImages.at(i))->GetData()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); directionImageContainer->InsertElement(directionImageContainer->Size(),itkImg); } catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } } // load reference directions ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); } catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } } // load/create mask image ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); if (maskImage.compare("")==0) { ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); itkMaskImage->SetSpacing( dirImg->GetSpacing() ); itkMaskImage->SetOrigin( dirImg->GetOrigin() ); itkMaskImage->SetDirection( dirImg->GetDirection() ); itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->Allocate(); itkMaskImage->FillBuffer(1); } else { mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::LoadDataNode(maskImage)->GetData()); mitk::CastToItkImage(mitkMaskImage, itkMaskImage); } // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); evaluationFilter->SetIgnoreMissingDirections(ignore); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); MITK_INFO << "writing " << outfilename; writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string logFile = outRoot; logFile.append("_ANGULAR_ERROR.csv"); ofstream file; file.open (logFile.c_str()); string sens = "Mean:"; sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(";\n"); sens.append("Median:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(";\n"); sens.append("Maximum:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(";\n"); sens.append("Minimum:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(";\n"); sens.append("STDEV:"); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; file.close(); MITK_INFO << "DONE"; } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } RegisterDiffusionMiniApp(PeaksAngularError); diff --git a/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp b/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp index be5bb44dea..b1994e085c 100644 --- a/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp +++ b/Modules/DiffusionImaging/MiniApps/QballReconstruction.cpp @@ -1,239 +1,244 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include "mitkBaseDataIOFactory.h" #include #include "mitkDiffusionImage.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include #include "ctkCommandLineParser.h" #include #include using namespace mitk; /** * Perform Q-ball reconstruction using a spherical harmonics basis */ int QballReconstruction(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("input", "i", ctkCommandLineParser::String, "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false); - parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output file", us::Any(), false); - parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "spherical harmonics order", 4, true); - parser.addArgument("b0Threshold", "t", ctkCommandLineParser::Int, "baseline image intensity threshold", 0, true); - parser.addArgument("lambda", "r", ctkCommandLineParser::Float, "ragularization factor lambda", 0.006, true); - parser.addArgument("csa", "csa", ctkCommandLineParser::Bool, "use constant solid angle consideration"); - parser.addArgument("outputCoeffs", "shc", ctkCommandLineParser::Bool, "output file containing the SH coefficients"); - parser.addArgument("mrtrix", "mb", ctkCommandLineParser::Bool, "use MRtrix compatible spherical harmonics definition"); + parser.addArgument("input", "i", ctkCommandLineParser::File, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false); + parser.addArgument("outFile", "o", ctkCommandLineParser::File, "Output file", "output file", us::Any(), false); + parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order", 4, true); + parser.addArgument("b0Threshold", "t", ctkCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0, true); + parser.addArgument("lambda", "r", ctkCommandLineParser::Float, "Lambda", "ragularization factor lambda", 0.006, true); + parser.addArgument("csa", "csa", ctkCommandLineParser::Bool, "Constant solid angle consideration", "use constant solid angle consideration"); + parser.addArgument("outputCoeffs", "shc", ctkCommandLineParser::Bool, "Output coefficients", "output file containing the SH coefficients"); + parser.addArgument("mrtrix", "mb", ctkCommandLineParser::Bool, "MRtrix", "use MRtrix compatible spherical harmonics definition"); + + parser.setCategory("Preprocessing Tools"); + parser.setTitle("Qball Reconstruction"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["input"]); std::string outfilename = us::any_cast(parsedArgs["outFile"]); outfilename = itksys::SystemTools::GetFilenamePath(outfilename)+"/"+itksys::SystemTools::GetFilenameWithoutExtension(outfilename); int threshold = 0; if (parsedArgs.count("b0Threshold")) threshold = us::any_cast(parsedArgs["b0Threshold"]); int shOrder = 4; if (parsedArgs.count("shOrder")) shOrder = us::any_cast(parsedArgs["shOrder"]); float lambda = 0.006; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); int normalization = 0; if (parsedArgs.count("csa") && us::any_cast(parsedArgs["csa"])) normalization = 6; bool outCoeffs = false; if (parsedArgs.count("outputCoeffs")) outCoeffs = us::any_cast(parsedArgs["outputCoeffs"]); bool mrTrix = false; if (parsedArgs.count("mrtrix")) mrTrix = us::any_cast(parsedArgs["mrtrix"]); try { MITK_INFO << "Loading image ..."; const std::string s1="", s2=""; std::vector infile = BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false ); DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); dwi->AverageRedundantGradients(0.001); mitk::QBallImage::Pointer image = mitk::QBallImage::New(); mitk::Image::Pointer coeffsImage = mitk::Image::New(); MITK_INFO << "SH order: " << shOrder; MITK_INFO << "lambda: " << lambda; MITK_INFO << "B0 threshold: " << threshold; switch ( shOrder ) { case 4: { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold( threshold ); filter->SetLambda(lambda); filter->SetUseMrtrixBasis(mrTrix); if (normalization==0) filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); else filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); break; } case 6: { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold( threshold ); filter->SetLambda(lambda); filter->SetUseMrtrixBasis(mrTrix); if (normalization==0) filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); else filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); break; } case 8: { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold( threshold ); filter->SetLambda(lambda); filter->SetUseMrtrixBasis(mrTrix); if (normalization==0) filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); else filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); break; } case 10: { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold( threshold ); filter->SetLambda(lambda); filter->SetUseMrtrixBasis(mrTrix); if (normalization==0) filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); else filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); break; } case 12: { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold( threshold ); filter->SetLambda(lambda); if (normalization==0) filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); else filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); break; } default: { MITK_INFO << "Supplied SH order not supported. Using default order of 4."; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold( threshold ); filter->SetLambda(lambda); filter->SetUseMrtrixBasis(mrTrix); if (normalization==0) filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); else filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); } } std::string coeffout = outfilename; coeffout += "_shcoeffs.nrrd"; outfilename += ".qbi"; MITK_INFO << "writing image " << outfilename; mitk::IOUtil::SaveBaseData(image, outfilename); if (outCoeffs) mitk::IOUtil::SaveImage(coeffsImage, coeffout); } catch ( itk::ExceptionObject &err) { MITK_INFO << "Exception: " << err; } catch ( std::exception err) { MITK_INFO << "Exception: " << err.what(); } catch ( ... ) { MITK_INFO << "Exception!"; } return EXIT_SUCCESS; } RegisterDiffusionMiniApp(QballReconstruction); diff --git a/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp b/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp index 4c3412402e..2605d62b90 100755 --- a/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp +++ b/Modules/DiffusionImaging/MiniApps/StreamlineTracking.cpp @@ -1,182 +1,187 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include #include #include #include #include #include #include #include "ctkCommandLineParser.h" #include #include int StreamlineTracking(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("input", "i", ctkCommandLineParser::StringList, "input tensor image (.dti)", us::Any(), false); - parser.addArgument("seed", "si", ctkCommandLineParser::String, "binary seed image", us::Any(), true); - parser.addArgument("mask", "mi", ctkCommandLineParser::String, "binary mask image", us::Any(), true); - parser.addArgument("faImage", "fai", ctkCommandLineParser::String, "FA image", us::Any(), true); - parser.addArgument("minFA", "fa", ctkCommandLineParser::Float, "minimum fractional anisotropy threshold", 0.15, true); - parser.addArgument("minCurv", "c", ctkCommandLineParser::Float, "minimum curvature radius in mm (default = 0.5*minimum-spacing)"); - parser.addArgument("stepSize", "s", ctkCommandLineParser::Float, "stepsize in mm (default = 0.1*minimum-spacing)"); - parser.addArgument("tendf", "f", ctkCommandLineParser::Float, "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0, true); - parser.addArgument("tendg", "g", ctkCommandLineParser::Float, "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0, true); - parser.addArgument("numSeeds", "n", ctkCommandLineParser::Int, "Number of seeds per voxel.", 1, true); - parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "minimum fiber length in mm", 20, true); - - parser.addArgument("interpolate", "ip", ctkCommandLineParser::Bool, "Use linear interpolation", false, true); - parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output fiber bundle (.fib)", us::Any(), false); + parser.addArgument("input", "i", ctkCommandLineParser::StringList, "Input image", "input tensor image (.dti)", us::Any(), false); + parser.addArgument("seed", "si", ctkCommandLineParser::File, "Seed image", "binary seed image", us::Any(), true); + parser.addArgument("mask", "mi", ctkCommandLineParser::File, "Mask", "binary mask image", us::Any(), true); + parser.addArgument("faImage", "fai", ctkCommandLineParser::File, "FA image", "FA image", us::Any(), true); + parser.addArgument("minFA", "fa", ctkCommandLineParser::Float, "Min. FA threshold", "minimum fractional anisotropy threshold", 0.15, true); + parser.addArgument("minCurv", "c", ctkCommandLineParser::Float, "Min. curvature radius", "minimum curvature radius in mm (default = 0.5*minimum-spacing)"); + parser.addArgument("stepSize", "s", ctkCommandLineParser::Float, "Step size", "step size in mm (default = 0.1*minimum-spacing)"); + parser.addArgument("tendf", "f", ctkCommandLineParser::Float, "Weight f", "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0, true); + parser.addArgument("tendg", "g", ctkCommandLineParser::Float, "Weight g", "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0, true); + parser.addArgument("numSeeds", "n", ctkCommandLineParser::Int, "Seeds per voxel", "Number of seeds per voxel.", 1, true); + parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "Min. fiber length", "minimum fiber length in mm", 20, true); + + parser.addArgument("interpolate", "ip", ctkCommandLineParser::Bool, "Interpolate", "Use linear interpolation", false, true); + parser.addArgument("outFile", "o", ctkCommandLineParser::String, "Output file", "output fiber bundle (.fib)", us::Any(), false); + + parser.setCategory("Fiber Tracking and Processing Methods"); + parser.setTitle("Streamline Tracking"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; ctkCommandLineParser::StringContainerType inputImages = us::any_cast(parsedArgs["input"]); string dtiFileName; string outFileName = us::any_cast(parsedArgs["outFile"]); float minFA = 0.15; float minCurv = -1; float stepSize = -1; float tendf = 1; float tendg = 0; float minLength = 20; int numSeeds = 1; bool interpolate = false; if (parsedArgs.count("minCurv")) minCurv = us::any_cast(parsedArgs["minCurv"]); if (parsedArgs.count("minFA")) minFA = us::any_cast(parsedArgs["minFA"]); if (parsedArgs.count("stepSize")) stepSize = us::any_cast(parsedArgs["stepSize"]); if (parsedArgs.count("tendf")) tendf = us::any_cast(parsedArgs["tendf"]); if (parsedArgs.count("tendg")) tendg = us::any_cast(parsedArgs["tendg"]); if (parsedArgs.count("minLength")) minLength = us::any_cast(parsedArgs["minLength"]); if (parsedArgs.count("numSeeds")) numSeeds = us::any_cast(parsedArgs["numSeeds"]); if (parsedArgs.count("interpolate")) interpolate = us::any_cast(parsedArgs["interpolate"]); try { typedef itk::StreamlineTrackingFilter< float > FilterType; FilterType::Pointer filter = FilterType::New(); mitk::Image::Pointer mitkImage = NULL; MITK_INFO << "Loading tensor images ..."; typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; dtiFileName = inputImages.at(0); for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(inputImages.at(i))->GetData()); mitk::TensorImage::Pointer img = dynamic_cast(mitk::IOUtil::LoadDataNode(inputImages.at(i))->GetData()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(img, itk_dti); filter->SetInput(i, itk_dti); } catch(...){ MITK_INFO << "could not load: " << inputImages.at(i); } } MITK_INFO << "Loading seed image ..."; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; mitk::Image::Pointer mitkSeedImage = NULL; if (parsedArgs.count("seed")) mitkSeedImage = mitk::IOUtil::LoadImage(us::any_cast(parsedArgs["seed"])); MITK_INFO << "Loading mask image ..."; mitk::Image::Pointer mitkMaskImage = NULL; if (parsedArgs.count("mask")) mitkMaskImage = mitk::IOUtil::LoadImage(us::any_cast(parsedArgs["mask"])); // instantiate tracker filter->SetSeedsPerVoxel(numSeeds); filter->SetFaThreshold(minFA); filter->SetMinCurvatureRadius(minCurv); filter->SetStepSize(stepSize); filter->SetF(tendf); filter->SetG(tendg); filter->SetInterpolate(interpolate); filter->SetMinTractLength(minLength); if (mitkSeedImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkSeedImage, mask); filter->SetSeedImage(mask); } if (mitkMaskImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMaskImage, mask); filter->SetMaskImage(mask); } filter->Update(); vtkSmartPointer fiberBundle = filter->GetFiberPolyData(); if ( fiberBundle->GetNumberOfLines()==0 ) { MITK_INFO << "No fibers reconstructed. Check parametrization."; return EXIT_FAILURE; } mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle); fib->SetReferenceImage(mitkImage); mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) { if ( (*it)->CanWriteBaseDataType(fib.GetPointer()) ) { (*it)->SetFileName( outFileName.c_str() ); (*it)->DoWrite( fib.GetPointer() ); } } } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_INFO << "DONE"; return EXIT_SUCCESS; } RegisterDiffusionMiniApp(StreamlineTracking); diff --git a/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp b/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp index debbad083a..77aec36383 100644 --- a/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp +++ b/Modules/DiffusionImaging/MiniApps/TensorDerivedMapsExtraction.cpp @@ -1,190 +1,195 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include #include "mitkImage.h" #include #include "mitkITKImageImport.h" #include #include #include #include #include "itkTensorDerivedMeasurementsFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "ctkCommandLineParser.h" #include #include #include typedef short DiffusionPixelType; typedef double TTensorPixelType; static void ExtractMapsAndSave(mitk::TensorImage::Pointer tensorImage, std::string filename, std::string postfix = "") { mitk::Image* image = dynamic_cast (tensorImage.GetPointer()); typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(image, itkvol); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); mitk::Image::Pointer map = mitk::Image::New(); // FA measurementsCalculator->SetMeasure(MeasurementsType::FA); measurementsCalculator->Update(); map->InitializeByItk( measurementsCalculator->GetOutput() ); map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() ); mitk::IOUtil::SaveImage(map, filename + "_dti_FA" + postfix + ".nrrd"); // MD measurementsCalculator->SetMeasure(MeasurementsType::MD); measurementsCalculator->Update(); map->InitializeByItk( measurementsCalculator->GetOutput() ); map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() ); mitk::IOUtil::SaveImage(map, filename + "_dti_MD" + postfix + ".nrrd"); // AD measurementsCalculator->SetMeasure(MeasurementsType::AD); measurementsCalculator->Update(); map->InitializeByItk( measurementsCalculator->GetOutput() ); map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() ); mitk::IOUtil::SaveImage(map, filename + "_dti_AD" + postfix + ".nrrd"); // CA measurementsCalculator->SetMeasure(MeasurementsType::CA); measurementsCalculator->Update(); map->InitializeByItk( measurementsCalculator->GetOutput() ); map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() ); mitk::IOUtil::SaveImage(map, filename + "_dti_CA" + postfix + ".nrrd"); // RA measurementsCalculator->SetMeasure(MeasurementsType::RA); measurementsCalculator->Update(); map->InitializeByItk( measurementsCalculator->GetOutput() ); map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() ); mitk::IOUtil::SaveImage(map, filename + "_dti_RA" + postfix + ".nrrd"); // RD measurementsCalculator->SetMeasure(MeasurementsType::RD); measurementsCalculator->Update(); map->InitializeByItk( measurementsCalculator->GetOutput() ); map->SetVolume( measurementsCalculator->GetOutput()->GetBufferPointer() ); mitk::IOUtil::SaveImage(map, filename + "_dti_RD" + postfix + ".nrrd"); } int TensorDerivedMapsExtraction(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("help", "h", ctkCommandLineParser::String, "Show this help text"); - parser.addArgument("input", "i", ctkCommandLineParser::String, "input dwi file", us::Any(),false); - parser.addArgument("out", "o", ctkCommandLineParser::String, "output folder and base name, e.g. /tmp/outPatient1 ", us::Any(),false); + parser.addArgument("help", "h", ctkCommandLineParser::String, "Help", "Show this help text"); + parser.addArgument("input", "i", ctkCommandLineParser::File, "Input file", "input dwi file", us::Any(),false); + parser.addArgument("out", "o", ctkCommandLineParser::String, "Output folder", "output folder and base name, e.g. /tmp/outPatient1 ", us::Any(),false); + + parser.setCategory("Diffusion Related Measures"); + parser.setTitle("Tensor Derived Maps Extraction"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << "\n\nMiniApp Description: \nPerforms tensor reconstruction on DWI file," << endl; std::cout << "and computes tensor derived measures." << endl; std::cout << "\n\n For out parameter /tmp/outPatient1 it will produce :"<< endl; std::cout << " /tmp/outPatient1_dti.dti , /tmp/outPatient1_dti_FA.nrrd, ..."<< endl; std::cout << "\n\n Parameters:"<< endl; std::cout << parser.helpText(); return EXIT_SUCCESS; } std::string inputFile = us::any_cast(parsedArgs["input"]); std::string baseFileName = us::any_cast(parsedArgs["out"]); std::string dtiFileName = "_dti.dti"; MITK_INFO << "BaseFileName: " << baseFileName; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputFile); mitk::DiffusionImage* diffusionImage = static_cast*>(inputImage.GetPointer()); if (diffusionImage == NULL) // does NULL pointer check make sense after static cast ? { MITK_ERROR << "Invalid Input Image. Must be DWI. Aborting."; return -1; } mitk::DiffusionImage* vols = dynamic_cast *> (inputImage.GetPointer()); typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientDirectionContainerType; GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = vols->GetDirections()->Begin(); it != vols->GetDirections()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() ); tensorReconstructionFilter->SetBValue(vols->GetReferenceBValue()); tensorReconstructionFilter->SetThreshold(50); tensorReconstructionFilter->Update(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer tensorImage = tensorReconstructionFilter->GetOutput(); tensorImage->SetDirection( vols->GetVectorImage()->GetDirection() ); mitk::TensorImage::Pointer tensorImageMitk = mitk::TensorImage::New(); tensorImageMitk->InitializeByItk(tensorImage.GetPointer()); tensorImageMitk->SetVolume( tensorImage->GetBufferPointer() ); itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); io->SetFileType( itk::ImageIOBase::Binary ); io->UseCompressionOn(); itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< double >, 3 > >::Pointer writer = itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< double >, 3 > >::New(); writer->SetInput(tensorReconstructionFilter->GetOutput()); writer->SetFileName(baseFileName + dtiFileName); writer->SetImageIO(io); writer->UseCompressionOn(); writer->Update(); ExtractMapsAndSave(tensorImageMitk,baseFileName); return EXIT_SUCCESS; } RegisterDiffusionMiniApp(TensorDerivedMapsExtraction); diff --git a/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp b/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp index d596931239..a9111bc8d5 100644 --- a/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp +++ b/Modules/DiffusionImaging/MiniApps/TensorReconstruction.cpp @@ -1,97 +1,102 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include "mitkBaseDataIOFactory.h" #include "mitkDiffusionImage.h" #include "mitkBaseData.h" #include #include #include #include #include "ctkCommandLineParser.h" #include using namespace mitk; /** * Convert files from one ending to the other */ int TensorReconstruction(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); - parser.addArgument("input", "i", ctkCommandLineParser::String, "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false); - parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output file", us::Any(), false); - parser.addArgument("b0Threshold", "t", ctkCommandLineParser::Int, "baseline image intensity threshold", 0, true); + parser.addArgument("input", "i", ctkCommandLineParser::File, "Input file", "input raw dwi (.dwi or .fsl/.fslgz)", us::Any(), false); + parser.addArgument("outFile", "o", ctkCommandLineParser::File, "Output file", "output file", us::Any(), false); + parser.addArgument("b0Threshold", "t", ctkCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0, true); + + parser.setCategory("Preprocessing Tools"); + parser.setTitle("Tensor Reconstruction"); + parser.setDescription(""); + parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["input"]); std::string outfilename = us::any_cast(parsedArgs["outFile"]); outfilename = itksys::SystemTools::GetFilenamePath(outfilename)+"/"+itksys::SystemTools::GetFilenameWithoutExtension(outfilename); outfilename += ".dti"; int threshold = 0; if (parsedArgs.count("b0Threshold")) threshold = us::any_cast(parsedArgs["b0Threshold"]); try { MITK_INFO << "Loading image ..."; const std::string s1="", s2=""; std::vector infile = BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false ); DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); MITK_INFO << "B0 threshold: " << threshold; typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, float > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetReferenceBValue()); filter->SetThreshold(threshold); filter->Update(); // Save tensor image MITK_INFO << "writing image " << outfilename; itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); io->SetFileType( itk::ImageIOBase::Binary ); io->UseCompressionOn(); itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< float >, 3 > >::Pointer writer = itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< float >, 3 > >::New(); writer->SetInput(filter->GetOutput()); writer->SetFileName(outfilename); writer->SetImageIO(io); writer->UseCompressionOn(); writer->Update(); } catch ( itk::ExceptionObject &err) { MITK_INFO << "Exception: " << err; } catch ( std::exception err) { MITK_INFO << "Exception: " << err.what(); } catch ( ... ) { MITK_INFO << "Exception!"; } return EXIT_SUCCESS; } RegisterDiffusionMiniApp(TensorReconstruction);