diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkCreation.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkCreation.cpp index bba63e81eb..d9646e0316 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkCreation.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkCreation.cpp @@ -1,119 +1,119 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // std includes #include // CTK includes #include "mitkCommandLineParser.h" // MITK includes #include "mitkConnectomicsNetworkCreator.h" #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Network Creation"); parser.setCategory("Connectomics"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "f", mitkCommandLineParser::String, "Input Tractogram", "input tractogram (.fib)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "p", mitkCommandLineParser::String, "Parcellation", "parcellation image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output network", "where to save the output (.cnf; .mat)", us::Any(), false); parser.addArgument("noCenterOfMass", "", mitkCommandLineParser::Bool, "No center of mass", "Do not use center of mass for node positions"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; //default values bool noCenterOfMass( false ); // parse command line arguments std::string fiberFilename = us::any_cast(parsedArgs["f"]); std::string parcellationFilename = us::any_cast(parsedArgs["p"]); std::string outputFilename = us::any_cast(parsedArgs["o"]); if (parsedArgs.count("noCenterOfMass")) noCenterOfMass = us::any_cast(parsedArgs["noCenterOfMass"]); try { // load fiber image std::vector fiberInfile = mitk::IOUtil::Load( fiberFilename); if( fiberInfile.empty() ) { std::string errorMessage = "Fiber Image at " + fiberFilename + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* fiberBaseData = fiberInfile.at(0); mitk::FiberBundle* fiberBundle = dynamic_cast( fiberBaseData ); // load parcellation std::vector parcellationInFile = mitk::IOUtil::Load( parcellationFilename); if( parcellationInFile.empty() ) { std::string errorMessage = "Parcellation at " + parcellationFilename + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* parcellationBaseData = parcellationInFile.at(0); mitk::Image* parcellationImage = dynamic_cast( parcellationBaseData ); // do creation mitk::ConnectomicsNetworkCreator::Pointer connectomicsNetworkCreator = mitk::ConnectomicsNetworkCreator::New(); connectomicsNetworkCreator->SetSegmentation( parcellationImage ); connectomicsNetworkCreator->SetFiberBundle( fiberBundle ); if( !noCenterOfMass ) { connectomicsNetworkCreator->CalculateCenterOfMass(); } connectomicsNetworkCreator->SetMappingStrategy(mitk::ConnectomicsNetworkCreator::MappingStrategy::EndElementPosition); connectomicsNetworkCreator->CreateNetworkFromFibersAndSegmentation(); mitk::ConnectomicsNetwork::Pointer network = connectomicsNetworkCreator->GetNetwork(); mitk::IOUtil::Save(network.GetPointer(), outputFilename ); return EXIT_SUCCESS; } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Connectomics/NetworkStatistics.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkStatistics.cpp index 411e993f70..758babe74d 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkStatistics.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkStatistics.cpp @@ -1,570 +1,570 @@ /*=================================================================== 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 #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Network Creation"); parser.setCategory("Connectomics"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input network", "input connectomics network (.cnf)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output file", "name of output file", us::Any(), false, false, false, mitkCommandLineParser::Output); 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",1); parser.addArgument("startDensity", "d", mitkCommandLineParser::Float, "Start Density", "Largest density for the range",1.0); parser.addArgument("thresholdStepSize", "t", mitkCommandLineParser::Int, "Step size threshold", "Distance of two adjacent thresholds",3); std::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["i"]); std::string outName = us::any_cast(parsedArgs["o"]); 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 { // load network std::vector networkFile = mitk::IOUtil::Load( networkName); if( networkFile.empty() ) { std::string errorMessage = "File at " + networkName + " could not be read. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } mitk::BaseData* networkBaseData = networkFile.at(0); mitk::ConnectomicsNetwork* network = dynamic_cast( networkBaseData ); if( !network ) { std::string errorMessage = "Read file at " + networkName + " could not be recognized as network. Aborting."; MITK_ERROR << errorMessage; return EXIT_FAILURE; } // streams std::stringstream globalHeaderStream; globalHeaderStream << "NumberOfVertices " << "NumberOfEdges " << "AverageDegree " << "ConnectionDensity " << "NumberOfConnectedComponents " << "AverageComponentSize " << "LargestComponentSize " << "RatioOfNodesInLargestComponent " << "HopPlotExponent " << "EffectiveHopDiameter " << "AverageClusteringCoefficientsC " << "AverageClusteringCoefficientsD " << "AverageClusteringCoefficientsE " << "AverageVertexBetweennessCentrality " << "AverageEdgeBetweennessCentrality " << "NumberOfIsolatedPoints " << "RatioOfIsolatedPoints " << "NumberOfEndPoints " << "RatioOfEndPoints " << "Diameter " << "Diameter90 " << "Radius " << "Radius90 " << "AverageEccentricity " << "AverageEccentricity90 " << "AveragePathLength " << "NumberOfCentralPoints " << "RatioOfCentralPoints " << "SpectralRadius " << "SecondLargestEigenValue " << "AdjacencyTrace " << "AdjacencyEnergy " << "LaplacianTrace " << "LaplacianEnergy " << "LaplacianSpectralGap " << "NormalizedLaplacianTrace " << "NormalizedLaplacianEnergy " << "NormalizedLaplacianNumberOf2s " << "NormalizedLaplacianNumberOf1s " << "NormalizedLaplacianNumberOf0s " << "NormalizedLaplacianLowerSlope " << "NormalizedLaplacianUpperSlope " << "SmallWorldness" << std::endl; std::stringstream localHeaderStream; std::stringstream regionalHeaderStream; std::stringstream globalDataStream; std::stringstream localDataStream; std::stringstream regionalDataStream; std::string globalOutName = outName + "_global.txt"; std::string localOutName = outName + "_local.txt"; std::string regionalOutName = outName + "_regional.txt"; bool firstRun( true ); // iterate over all three possible methods for(unsigned int method( 0 ); method < 3; method++) { // 0 - Random removal threshold // 1 - Largest density below threshold // 2 - Threshold based // iterate over possible targets for( int step = 0; step < granularity; ++step ) { double targetValue( 0.0 ); 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; mitk::LocaleSwitch localeSwitch("C"); 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(std::size_t 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( std::size_t 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; // count number of connections and fibers between regions std::map< std::string, std::vector >::iterator loopRegionsIterator; for (loopRegionsIterator = parsedRegionsIterator; loopRegionsIterator != parsedRegions.end(); loopRegionsIterator++) { int numberConnections(0), possibleConnections(0); double summedFiberCount(0.0); std::vector loopLabelsVector = loopRegionsIterator->second; std::string loopName = loopRegionsIterator->first; for (std::size_t loop(0); loop < regionLabelsVector.size(); loop++) { if (thresholdedNetwork->CheckForLabel(regionLabelsVector.at(loop))) { for (std::size_t innerLoop(0); innerLoop < loopLabelsVector.size(); innerLoop++) { if (thresholdedNetwork->CheckForLabel(loopLabelsVector.at(loop))) { bool exists = thresholdedNetwork->EdgeExists( labelToIdMap.find(regionLabelsVector.at(loop))->second, labelToIdMap.find(loopLabelsVector.at(innerLoop))->second); possibleConnections++; if (exists) { numberConnections++; summedFiberCount += thresholdedNetwork->GetEdge( labelToIdMap.find(regionLabelsVector.at(loop))->second, labelToIdMap.find(loopLabelsVector.at(innerLoop))->second).fiber_count; } } else { MITK_ERROR << "Illegal label. Label: \"" << loopLabelsVector.at(loop) << "\" not found."; } } } else { MITK_ERROR << "Illegal label. Label: \"" << regionLabelsVector.at(loop) << "\" not found."; } } if (firstRun) { regionalHeaderStream << " " << regionName << "_" << loopName << "_Connections " << " " << regionName << "_" << loopName << "_possibleConnections " << " " << regionName << "_" << loopName << "_ConnectingFibers"; } regionalDataStream << " " << numberConnections << " " << possibleConnections << " " << summedFiberCount; } } } 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 (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FiberClustering.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberClustering.cpp index 94ca3fb900..0ddcd5aff6 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberClustering.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberClustering.cpp @@ -1,273 +1,273 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Clustering"); parser.setCategory("Fiber Processing"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle (.fib; .trk; .tck)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("cluster_size", "", mitkCommandLineParser::Int, "Cluster size:", "", 10); parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 12); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. fibers per cluster:", "", 1); parser.addArgument("max_clusters", "", mitkCommandLineParser::Int, "Max. clusters:", ""); parser.addArgument("merge_clusters", "", mitkCommandLineParser::Float, "Merge clusters:", "Set to 0 to avoid merging and to -1 to use the original cluster size", -1.0); parser.addArgument("output_centroids", "", mitkCommandLineParser::Bool, "Output centroids:", ""); parser.addArgument("only_centroids", "", mitkCommandLineParser::Bool, "Output only centroids:", ""); parser.addArgument("merge_centroids", "", mitkCommandLineParser::Bool, "Merge centroids:", ""); parser.addArgument("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN; EU_STD; EU_MAX; ANAT; MAP; LENGTH", std::string("EU_MEAN")); parser.addArgument("metric_weights", "", mitkCommandLineParser::StringList, "Metric weights:", "add one float weight for each used metric"); parser.addArgument("input_centroids", "", mitkCommandLineParser::String, "Input centroids:", "", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("scalar_map", "", mitkCommandLineParser::String, "Scalar map:", "", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("parcellation", "", mitkCommandLineParser::String, "Parcellation:", "", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("file_ending", "", mitkCommandLineParser::String, "File ending:", ""); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); bool only_centroids = false; if (parsedArgs.count("only_centroids")) only_centroids = us::any_cast(parsedArgs["only_centroids"]); bool merge_centroids = false; if (parsedArgs.count("merge_centroids")) merge_centroids = us::any_cast(parsedArgs["merge_centroids"]); int cluster_size = 10; if (parsedArgs.count("cluster_size")) cluster_size = us::any_cast(parsedArgs["cluster_size"]); int fiber_points = 12; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); int min_fibers = 1; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); int max_clusters = 0; if (parsedArgs.count("max_clusters")) max_clusters = us::any_cast(parsedArgs["max_clusters"]); float merge_clusters = -1.0; if (parsedArgs.count("merge_clusters")) merge_clusters = us::any_cast(parsedArgs["merge_clusters"]); bool output_centroids = false; if (parsedArgs.count("output_centroids")) output_centroids = us::any_cast(parsedArgs["output_centroids"]); std::vector< std::string > metric_strings = {"EU_MEAN"}; if (parsedArgs.count("metrics")) metric_strings = us::any_cast(parsedArgs["metrics"]); std::vector< std::string > metric_weights = {"1.0"}; if (parsedArgs.count("metric_weights")) metric_weights = us::any_cast(parsedArgs["metric_weights"]); std::string input_centroids = ""; if (parsedArgs.count("input_centroids")) input_centroids = us::any_cast(parsedArgs["input_centroids"]); std::string scalar_map = ""; if (parsedArgs.count("scalar_map")) scalar_map = us::any_cast(parsedArgs["scalar_map"]); std::string parcellation = ""; if (parsedArgs.count("parcellation")) parcellation = us::any_cast(parsedArgs["parcellation"]); std::string file_ending = ".fib"; if (parsedArgs.count("file_ending")) file_ending = us::any_cast(parsedArgs["file_ending"]); if (metric_strings.size()!=metric_weights.size()) { MITK_INFO << "Each metric needs an associated metric weight!"; return EXIT_FAILURE; } try { typedef itk::Image< float, 3 > FloatImageType; typedef itk::Image< short, 3 > ShortImageType; mitk::FiberBundle::Pointer fib = LoadFib(inFileName); float max_d = 0; int i=1; std::vector< float > distances; while (max_d < fib->GetGeometry()->GetDiagonalLength()/2) { distances.push_back(cluster_size*i); max_d = cluster_size*i; ++i; } itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances(distances); clusterer->SetTractogram(fib); if (input_centroids!="") { mitk::FiberBundle::Pointer in_centroids = LoadFib(input_centroids); clusterer->SetInCentroids(in_centroids); } std::vector< mitk::ClusteringMetric* > metrics; int mc = 0; for (auto m : metric_strings) { float w = boost::lexical_cast(metric_weights.at(mc)); MITK_INFO << "Metric: " << m << " (w=" << w << ")"; if (m=="EU_MEAN") metrics.push_back({new mitk::ClusteringMetricEuclideanMean()}); else if (m=="EU_STD") metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); else if (m=="EU_MAX") metrics.push_back({new mitk::ClusteringMetricEuclideanMax()}); else if (m=="ANGLES") metrics.push_back({new mitk::ClusteringMetricInnerAngles()}); else if (m=="LENGTH") metrics.push_back({new mitk::ClusteringMetricLength()}); else if (m=="MAP" && scalar_map!="") { mitk::Image::Pointer mitk_map = mitk::IOUtil::Load(scalar_map); if (mitk_map->GetDimension()==3) { FloatImageType::Pointer itk_map = FloatImageType::New(); mitk::CastToItkImage(mitk_map, itk_map); mitk::ClusteringMetricScalarMap* metric = new mitk::ClusteringMetricScalarMap(); metric->SetImages({itk_map}); metric->SetScale(distances.at(0)); metrics.push_back(metric); } } else if (m=="ANAT" && parcellation!="") { mitk::Image::Pointer mitk_map = mitk::IOUtil::Load(parcellation); if (mitk_map->GetDimension()==3) { ShortImageType::Pointer itk_map = ShortImageType::New(); mitk::CastToItkImage(mitk_map, itk_map); mitk::ClusteringMetricAnatomic* metric = new mitk::ClusteringMetricAnatomic(); metric->SetParcellations({itk_map}); metrics.push_back(metric); } } metrics.back()->SetScale(w); mc++; } if (metrics.empty()) { MITK_INFO << "No metric selected!"; return EXIT_FAILURE; } clusterer->SetMetrics(metrics); clusterer->SetMergeDuplicateThreshold(merge_clusters); clusterer->SetNumPoints(fiber_points); clusterer->SetMaxClusters(max_clusters); clusterer->SetMinClusterSize(min_fibers); clusterer->Update(); std::vector tracts = clusterer->GetOutTractograms(); std::vector centroids = clusterer->GetOutCentroids(); MITK_INFO << "Saving clusters"; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect if (!only_centroids) for (unsigned int i=0; i(i) + file_ending); if (output_centroids && !merge_centroids) { for (unsigned int i=0; i(i) + file_ending); } else if (output_centroids) { mitk::FiberBundle::Pointer centroid = mitk::FiberBundle::New(); centroid = centroid->AddBundles(centroids); mitk::IOUtil::Save(centroid, out_root + ist::GetFilenameWithoutExtension(inFileName) + "_Centroids" + file_ending); } std::cout.rdbuf (old); // <-- restore } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FiberDirectionExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberDirectionExtraction.cpp index 0d3fbe4f2a..4bb955b3af 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberDirectionExtraction.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberDirectionExtraction.cpp @@ -1,183 +1,183 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include /*! \brief Extract principal fiber directions from a tractogram */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Direction Extraction"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Extract principal fiber directions from a tractogram"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("mask", "", mitkCommandLineParser::String, "Mask:", "mask image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("athresh", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); parser.addArgument("peakthresh", "", mitkCommandLineParser::Float, "Peak size threshold:", "peak size threshold relative to largest peak in voxel", 0.2, true); parser.addArgument("only_mask_geometry", "", mitkCommandLineParser::Bool, "Only mask geometry:", "don't use content of mask image, only use it's geometry", false); parser.addArgument("verbose", "", mitkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results"); parser.addArgument("numdirs", "", mitkCommandLineParser::Int, "Max. num. directions:", "maximum number of fibers per voxel", 3, true); parser.addArgument("normalization", "", mitkCommandLineParser::Int, "Normalization method:", "1=global maximum; 2=single vector; 3=voxel-wise maximum", 1); parser.addArgument("file_ending", "", mitkCommandLineParser::String, "Image type:", ".nrrd; .nii; .nii.gz", std::string(".nii.gz")); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string fibFile = us::any_cast(parsedArgs["i"]); std::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"]); std::string outRoot = us::any_cast(parsedArgs["o"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); bool only_mask_geometry = false; if (parsedArgs.count("only_mask_geometry")) only_mask_geometry = us::any_cast(parsedArgs["only_mask_geometry"]); int maxNumDirs = 3; if (parsedArgs.count("numdirs")) maxNumDirs = us::any_cast(parsedArgs["numdirs"]); int normalization = 1; if (parsedArgs.count("normalization")) normalization = us::any_cast(parsedArgs["normalization"]); std::string file_ending = ".nii.gz"; if (parsedArgs.count("file_ending")) file_ending = us::any_cast(parsedArgs["file_ending"]); try { typedef itk::Image ItkUcharImgType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(fibFile); // load/create mask image ItkUcharImgType::Pointer itkMaskImage = nullptr; if (maskImage.compare("")!=0) { std::cout << "Using mask image"; itkMaskImage = ItkUcharImgType::New(); mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::Load(maskImage); mitk::CastToItkImage(mitkMaskImage, itkMaskImage); } // extract directions from fiber bundle itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetOnlyUseMaskGeometry(only_mask_geometry); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*itk::Math::pi/180)); switch (normalization) { case 1: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::GLOBAL_MAX); break; case 2: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::SINGLE_VEC_NORM); break; case 3: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::MAX_VEC_NORM); break; } fOdfFilter->SetSizeThreshold(peakThreshold); fOdfFilter->SetMaxNumDirections(maxNumDirs); fOdfFilter->Update(); mitk::LocaleSwitch localeSwitch("C"); { itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = fOdfFilter->GetDirectionImage(); typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; WriterType::Pointer writer = WriterType::New(); std::string outfilename = outRoot; outfilename.append("_DIRECTIONS"); outfilename.append(file_ending); writer->SetFileName(outfilename.c_str()); writer->SetInput(itkImg); writer->Update(); } if (verbose) { // write num direction image ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); std::string outfilename = outRoot; outfilename.append("_NUM_DIRECTIONS"); outfilename.append(file_ending); writer->SetFileName(outfilename.c_str()); writer->SetInput(numDirImage); writer->Update(); } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FiberExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtraction.cpp index 16a0c2fb19..9eb0dbeedc 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtraction.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtraction.cpp @@ -1,151 +1,151 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include /*! \brief Extract fibers from a tractogram using planar figure ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Extraction"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Extract fibers from a tractogram using planar figure ROIs"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); parser.addArgument("planfirgure1", "pf1", mitkCommandLineParser::String, "Figure 1:", "first planar figure ROI", us::Any(), false); parser.addArgument("planfirgure2", "pf2", mitkCommandLineParser::String, "Figure 2:", "second planar figure ROI", us::Any()); parser.addArgument("operation", "op", mitkCommandLineParser::String, "Operation:", "logical operation (AND; OR; NOT)", us::Any()); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["i"]); std::string outFib = us::any_cast(parsedArgs["o"]); std::string pf1_path = us::any_cast(parsedArgs["planfirgure1"]); std::string operation(""); std::string pf2_path(""); if (parsedArgs.count("operation")) { operation = us::any_cast(parsedArgs["operation"]); if (parsedArgs.count("planfirgure2") && (operation=="AND" || operation=="OR")) pf2_path = us::any_cast(parsedArgs["planfirgure2"]); } try { // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(inFib); mitk::FiberBundle::Pointer result; mitk::StandaloneDataStorage::Pointer storage = mitk::StandaloneDataStorage::New(); auto data = mitk::IOUtil::Load(pf1_path)[0]; auto input1 = mitk::DataNode::New(); input1->SetData(data); if (input1.IsNotNull()) { mitk::PlanarFigureComposite::Pointer pfc = mitk::PlanarFigureComposite::New(); mitk::DataNode::Pointer pfcNode = mitk::DataNode::New(); pfcNode->SetData(pfc); mitk::DataStorage::SetOfObjects::Pointer set1 = mitk::DataStorage::SetOfObjects::New(); set1->push_back(pfcNode); storage->Add(pfcNode); auto input2 = mitk::DataNode::New(); if (!pf2_path.empty()) { data = mitk::IOUtil::Load(pf2_path)[0]; input2->SetData(data); } if (operation.empty()) { result = inputTractogram->ExtractFiberSubset(input1, nullptr); } else if (operation=="NOT") { pfc->setOperationType(mitk::PlanarFigureComposite::NOT); storage->Add(input1, set1); result = inputTractogram->ExtractFiberSubset(pfcNode, storage); } else if (operation=="AND" && input2.IsNotNull()) { pfc->setOperationType(mitk::PlanarFigureComposite::AND); storage->Add(input1, set1); storage->Add(input2, set1); result = inputTractogram->ExtractFiberSubset(pfcNode, storage); } else if (operation=="OR" && input2.IsNotNull()) { pfc->setOperationType(mitk::PlanarFigureComposite::OR); storage->Add(input1, set1); storage->Add(input2, set1); result = inputTractogram->ExtractFiberSubset(pfcNode, storage); } else { std::cout << "Could not process input:"; std::cout << pf1_path; std::cout << pf2_path; std::cout << operation; } } if (result.IsNotNull()) mitk::IOUtil::Save(result, outFib); else std::cout << "No valid fiber bundle extracted."; } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp index 5c55d70743..a464985d16 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp @@ -1,276 +1,276 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = mitk::IOUtil::Load(filename); ItkFloatImgType::Pointer itk_image = ItkFloatImgType::New(); mitk::CastToItkImage(img, itk_image); return itk_image; } /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Extraction With ROI Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Extract fibers from a tractogram using binary image ROIs"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false); parser.addArgument("rois", "", mitkCommandLineParser::StringList, "ROI images:", "ROI images", us::Any(), false); parser.endGroup(); parser.beginGroup("2. Label based extraction:"); parser.addArgument("split_labels", "", mitkCommandLineParser::Bool, "Split labels:", "output a separate tractogram for each label-->label tract", false); parser.addArgument("skip_self_connections", "", mitkCommandLineParser::Bool, "Skip self connections:", "ignore streamlines between two identical labels", false); parser.addArgument("all_labels", "", mitkCommandLineParser::Bool, "All labels:", "use all labels (0 is excluded)", false); parser.addArgument("labels", "", mitkCommandLineParser::StringList, "Labels:", "positive means roi image value in labels vector", us::Any()); parser.addArgument("start_labels", "", mitkCommandLineParser::StringList, "Start Labels:", "use separate start and end labels instead of one mixed set", us::Any()); parser.addArgument("end_labels", "", mitkCommandLineParser::StringList, "End Labels:", "use separate start and end labels instead of one mixed set", us::Any()); parser.addArgument("paired", "", mitkCommandLineParser::Bool, "Paired:", "start and end label list are paired", false); parser.endGroup(); parser.beginGroup("3. Misc:"); parser.addArgument("both_ends", "", mitkCommandLineParser::Bool, "Both ends:", "Fibers are extracted if both endpoints are located in the ROI.", false); parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "Extract by overlap, not by endpoints. Extract fibers that overlap to at least the provided factor (0-1) with the ROI.", -1); parser.addArgument("invert", "", mitkCommandLineParser::Bool, "Invert:", "get streamlines not positive for any of the ROI images", false); parser.addArgument("output_negatives", "", mitkCommandLineParser::Bool, "Negatives:", "output negatives", false); parser.addArgument("interpolate", "", mitkCommandLineParser::Bool, "Interpolate:", "interpolate ROI images (only for endpoint based extraction)", false); parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", 0.5); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. num. fibers:", "discard positive tracts with less fibers", 0); parser.addArgument("split_rois", "", mitkCommandLineParser::Bool, "Split ROIs:", "output a separate tractogram for each ROI", false); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["i"]); std::string outFib = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType roi_files = us::any_cast(parsedArgs["rois"]); bool both_ends = false; if (parsedArgs.count("both_ends")) both_ends = us::any_cast(parsedArgs["both_ends"]); bool invert = false; if (parsedArgs.count("invert")) invert = us::any_cast(parsedArgs["invert"]); unsigned int min_fibers = 0; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); bool all_labels = false; if (parsedArgs.count("all_labels")) all_labels = us::any_cast(parsedArgs["all_labels"]); bool split_labels = false; if (parsedArgs.count("split_labels")) split_labels = us::any_cast(parsedArgs["split_labels"]); bool split_rois = false; if (parsedArgs.count("split_rois")) split_rois = us::any_cast(parsedArgs["split_rois"]); bool skip_self_connections = false; if (parsedArgs.count("skip_self_connections")) skip_self_connections = us::any_cast(parsedArgs["skip_self_connections"]); bool output_negatives = false; if (parsedArgs.count("output_negatives")) output_negatives = us::any_cast(parsedArgs["output_negatives"]); float overlap_fraction = -1; if (parsedArgs.count("overlap_fraction")) overlap_fraction = us::any_cast(parsedArgs["overlap_fraction"]); bool any_point = false; if (overlap_fraction>=0) any_point = true; bool interpolate = false; if (parsedArgs.count("interpolate")) interpolate = us::any_cast(parsedArgs["interpolate"]); float threshold = 0.5; if (parsedArgs.count("threshold")) threshold = us::any_cast(parsedArgs["threshold"]); mitkCommandLineParser::StringContainerType labels; if (parsedArgs.count("labels")) labels = us::any_cast(parsedArgs["labels"]); mitkCommandLineParser::StringContainerType start_labels; if (parsedArgs.count("start_labels")) start_labels = us::any_cast(parsedArgs["start_labels"]); mitkCommandLineParser::StringContainerType end_labels; if (parsedArgs.count("end_labels")) end_labels = us::any_cast(parsedArgs["end_labels"]); bool paired = false; if (parsedArgs.count("paired")) paired = us::any_cast(parsedArgs["paired"]); try { // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(inFib); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect std::vector< ItkFloatImgType::Pointer > roi_images; std::vector< std::string > roi_names; for (std::size_t i=0; i roi_images2; for (auto roi : roi_images) roi_images2.push_back(roi); std::vector< unsigned short > short_labels; for (auto l : labels) short_labels.push_back(boost::lexical_cast(l)); std::vector< unsigned short > short_start_labels; for (auto l : start_labels) short_start_labels.push_back(boost::lexical_cast(l)); std::vector< unsigned short > short_end_labels; for (auto l : end_labels) short_end_labels.push_back(boost::lexical_cast(l)); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(inputTractogram); extractor->SetRoiImages(roi_images2); extractor->SetRoiImageNames(roi_names); extractor->SetOverlapFraction(overlap_fraction); extractor->SetBothEnds(both_ends); extractor->SetInterpolate(interpolate); extractor->SetThreshold(threshold); extractor->SetLabels(short_labels); extractor->SetStartLabels(short_start_labels); extractor->SetEndLabels(short_end_labels); extractor->SetSplitLabels(split_labels); extractor->SetSplitByRoi(split_rois); extractor->SetMinFibersPerTract(min_fibers); extractor->SetSkipSelfConnections(skip_self_connections); extractor->SetPairedStartEndLabels(paired); if (!any_point) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); if (all_labels || short_labels.size()>0 || short_start_labels.size()>0 || short_end_labels.size()>0) extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); extractor->Update(); std::string ext = itksys::SystemTools::GetFilenameExtension(outFib); if (ext.empty()) ext = ".trk"; outFib = itksys::SystemTools::GetFilenamePath(outFib) + '/' + itksys::SystemTools::GetFilenameWithoutExtension(outFib); if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib + ext); else { int c = 0; std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); for (auto fib : extractor->GetPositives()) { std::string l = positive_labels.at(c); if (l.size()>0) mitk::IOUtil::Save(fib, outFib + "_" + l + ext); else mitk::IOUtil::Save(fib, outFib + ext); ++c; } } if (output_negatives) { invert = !invert; if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib + "_negatives" + ext); else { int c = 0; std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); for (auto fib : extractor->GetPositives()) { std::string l = positive_labels.at(c); if (l.size()>0) mitk::IOUtil::Save(fib, outFib + "_" + l + "_negatives" + ext); else mitk::IOUtil::Save(fib, outFib + "_negatives" + ext); ++c; } } } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FiberJoin.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberJoin.cpp index 147bb3db9d..fb9bcdd533 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberJoin.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberJoin.cpp @@ -1,100 +1,100 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Join multiple tractograms */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Join"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Join multiple tractograms"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::StringList, "Input:", "input tractograms", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType inFibs = us::any_cast(parsedArgs["i"]); std::string outFib = us::any_cast(parsedArgs["o"]); if (inFibs.size()<=1) { std::cout << "More than one input tractogram required!"; return EXIT_FAILURE; } try { std::vector< mitk::FiberBundle::Pointer > tractograms; mitk::FiberBundle::Pointer result = LoadFib(inFibs.at(0)); for (std::size_t i=1; iAddBundles(tractograms); mitk::IOUtil::Save(result, outFib); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FiberProcessing.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberProcessing.cpp index b9ea3f9d5a..648fe2042e 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberProcessing.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberProcessing.cpp @@ -1,301 +1,301 @@ /*=================================================================== 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 #include #include #include mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Modify input tractogram: fiber resampling, compression, pruning and transformation. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Processing"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Modify input tractogram: fiber resampling, compression, pruning and transformation."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "Input fiber bundle (.fib, .trk, .tck)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "Output fiber bundle (.fib, .trk)", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.endGroup(); parser.beginGroup("2. Resampling:"); parser.addArgument("spline_resampling", "", mitkCommandLineParser::Float, "Spline resampling:", "Resample fiber using splines with the given point distance (in mm)"); parser.addArgument("linear_resampling", "", mitkCommandLineParser::Float, "Linear resampling:", "Resample fiber linearly with the given point distance (in mm)"); parser.addArgument("num_resampling", "", mitkCommandLineParser::Int, "Num. fiber points resampling:", "Resample all fibers to the given number of points"); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "Compress fiber using the given error threshold (in mm)"); parser.endGroup(); parser.beginGroup("3. Filtering:"); parser.addArgument("min_length", "", mitkCommandLineParser::Float, "Minimum length:", "Minimum fiber length (in mm)"); parser.addArgument("max_length", "", mitkCommandLineParser::Float, "Maximum length:", "Maximum fiber length (in mm)"); parser.addArgument("max_angle", "", mitkCommandLineParser::Float, "Maximum angle:", "Maximum angular STDEV (in degree) over given distance"); parser.addArgument("max_angle_dist", "", mitkCommandLineParser::Float, "Distance:", "Distance in mm", 10); parser.addArgument("remove", "", mitkCommandLineParser::Bool, "Remove fibers exceeding curvature threshold:", "If false, only the high curvature parts are removed"); parser.addArgument("subsample", "", mitkCommandLineParser::Float, "Randomly select fraction of streamlines:", "Randomly select the specified fraction of streamlines from the input tractogram"); parser.addArgument("random_subsample", "", mitkCommandLineParser::Bool, "Randomly seed subsampling:", "Randomly seed subsampling. Else, use seed 0."); parser.endGroup(); parser.beginGroup("4. Transformation:"); parser.addArgument("mirror", "", mitkCommandLineParser::Int, "Invert coordinates:", "Invert fiber coordinates XYZ (e.g. 010 to invert y-coordinate of each fiber point)"); parser.addArgument("rotate_x", "", mitkCommandLineParser::Float, "Rotate x-axis:", "Rotate around x-axis (in deg)"); parser.addArgument("rotate_y", "", mitkCommandLineParser::Float, "Rotate y-axis:", "Rotate around y-axis (in deg)"); parser.addArgument("rotate_z", "", mitkCommandLineParser::Float, "Rotate z-axis:", "Rotate around z-axis (in deg)"); parser.addArgument("scale_x", "", mitkCommandLineParser::Float, "Scale x-axis:", "Scale in direction of x-axis"); parser.addArgument("scale_y", "", mitkCommandLineParser::Float, "Scale y-axis:", "Scale in direction of y-axis"); parser.addArgument("scale_z", "", mitkCommandLineParser::Float, "Scale z-axis", "Scale in direction of z-axis"); parser.addArgument("translate_x", "", mitkCommandLineParser::Float, "Translate x-axis:", "Translate in direction of x-axis (in mm)"); parser.addArgument("translate_y", "", mitkCommandLineParser::Float, "Translate y-axis:", "Translate in direction of y-axis (in mm)"); parser.addArgument("translate_z", "", mitkCommandLineParser::Float, "Translate z-axis:", "Translate in direction of z-axis (in mm)"); parser.endGroup(); parser.beginGroup("5. Remove fiber parts:"); parser.addArgument("remove_inside", "", mitkCommandLineParser::Bool, "Remove fibers inside mask:", "remove fibers inside mask"); parser.addArgument("remove_outside", "", mitkCommandLineParser::Bool, "Remove fibers outside mask:", "remove fibers outside mask"); parser.addArgument("mask", "", mitkCommandLineParser::String, "Mask image:", "mask image", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; bool remove_outside = false; if (parsedArgs.count("remove_outside")) remove_outside = us::any_cast(parsedArgs["remove_outside"]); bool remove_inside = false; if (!remove_outside && parsedArgs.count("remove_inside")) remove_inside = us::any_cast(parsedArgs["remove_inside"]); typedef itk::Image< unsigned char, 3 > UcharImageType; UcharImageType::Pointer mask = nullptr; if (remove_inside || remove_outside) { if (parsedArgs.count("mask")) mask = mitk::DiffusionDataIOHelper::load_itk_image< UcharImageType >(us::any_cast(parsedArgs["mask"])); else { MITK_INFO << "Mask needed to remove fibers inside or outside mask!"; return EXIT_FAILURE; } } bool remove = false; if (parsedArgs.count("remove")) remove = us::any_cast(parsedArgs["remove"]); bool random_subsample = false; if (parsedArgs.count("random_subsample")) random_subsample = us::any_cast(parsedArgs["random_subsample"]); float spline_resampling = -1; if (parsedArgs.count("spline_resampling")) spline_resampling = us::any_cast(parsedArgs["spline_resampling"]); float linear_resampling = -1; if (parsedArgs.count("linear_resampling")) linear_resampling = us::any_cast(parsedArgs["linear_resampling"]); int num_resampling = -1; if (parsedArgs.count("num_resampling")) num_resampling = us::any_cast(parsedArgs["num_resampling"]); float subsample = -1; if (parsedArgs.count("subsample")) subsample = us::any_cast(parsedArgs["subsample"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); float minFiberLength = -1; if (parsedArgs.count("min_length")) minFiberLength = us::any_cast(parsedArgs["min_length"]); float maxFiberLength = -1; if (parsedArgs.count("max_length")) maxFiberLength = us::any_cast(parsedArgs["max_length"]); float max_angle_dist = 10; if (parsedArgs.count("max_angle_dist")) max_angle_dist = us::any_cast(parsedArgs["max_angle_dist"]); float maxAngularDev = -1; if (parsedArgs.count("max_angle")) maxAngularDev = us::any_cast(parsedArgs["max_angle"]); int axis = 0; if (parsedArgs.count("mirror")) axis = us::any_cast(parsedArgs["mirror"]); float rotateX = 0; if (parsedArgs.count("rotate_x")) rotateX = us::any_cast(parsedArgs["rotate_x"]); float rotateY = 0; if (parsedArgs.count("rotate_y")) rotateY = us::any_cast(parsedArgs["rotate_y"]); float rotateZ = 0; if (parsedArgs.count("rotate_z")) rotateZ = us::any_cast(parsedArgs["rotate_z"]); float scaleX = 0; if (parsedArgs.count("scale_x")) scaleX = us::any_cast(parsedArgs["scale_x"]); float scaleY = 0; if (parsedArgs.count("scale_y")) scaleY = us::any_cast(parsedArgs["scale_y"]); float scaleZ = 0; if (parsedArgs.count("scale_z")) scaleZ = us::any_cast(parsedArgs["scale_z"]); float translateX = 0; if (parsedArgs.count("translate_x")) translateX = us::any_cast(parsedArgs["translate_x"]); float translateY = 0; if (parsedArgs.count("translate_y")) translateY = us::any_cast(parsedArgs["translate_y"]); float translateZ = 0; if (parsedArgs.count("translate_z")) translateZ = us::any_cast(parsedArgs["translate_z"]); std::string inFileName = us::any_cast(parsedArgs["i"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { mitk::FiberBundle::Pointer fib = LoadFib(inFileName); if (subsample>0) fib = fib->SubsampleFibers(subsample, random_subsample); if (maxAngularDev>0) { auto filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(maxAngularDev); filter->SetDistance(max_angle_dist); filter->SetRemoveFibers(remove); filter->Update(); fib = filter->GetOutputFiberBundle(); } if (minFiberLength>0) fib->RemoveShortFibers(minFiberLength); if (maxFiberLength>0) fib->RemoveLongFibers(maxFiberLength); if (spline_resampling>0) fib->ResampleSpline(spline_resampling); if (linear_resampling>0) fib->ResampleLinear(linear_resampling); if (num_resampling>0) fib->ResampleToNumPoints(num_resampling); if (compress>0) fib->Compress(compress); if ( mask.IsNotNull() ) { if (remove_outside) fib = fib->RemoveFibersOutside(mask, false); else if (remove_inside) fib = fib->RemoveFibersOutside(mask, true); } if (axis/100==1) fib->MirrorFibers(0); if ((axis%100)/10==1) fib->MirrorFibers(1); if (axis%10==1) fib->MirrorFibers(2); if (rotateX > 0 || rotateY > 0 || rotateZ > 0){ std::cout << "Rotate " << rotateX << " " << rotateY << " " << rotateZ; fib->RotateAroundAxis(rotateX, rotateY, rotateZ); } if (translateX > 0 || translateY > 0 || translateZ > 0){ fib->TranslateFibers(translateX, translateY, translateZ); } if (scaleX > 0 || scaleY > 0 || scaleZ > 0) fib->ScaleFibers(scaleX, scaleY, scaleZ); mitk::IOUtil::Save(fib.GetPointer(), outFileName ); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp index fba39ecc48..bd9b319e7d 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp @@ -1,321 +1,321 @@ /*=================================================================== 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 typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; /*! \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092). */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fit Fibers To Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Assigns a weight to each fiber in order to optimally explain the input peak image"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkCommandLineParser::StringList, "Input tractograms:", "input tractograms (files or folder)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "i2", mitkCommandLineParser::String, "Input image:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("max_iter", "", mitkCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20); parser.addArgument("bundle_based", "", mitkCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false); parser.addArgument("min_g", "", mitkCommandLineParser::Float, "Min. g:", "lower termination threshold for gradient magnitude", 1e-5); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 1.0); parser.addArgument("save_res", "", mitkCommandLineParser::Bool, "Save Residuals:", "save residual images", false); parser.addArgument("save_weights", "", mitkCommandLineParser::Bool, "Save Weights:", "save fiber weights in a separate text file", false); parser.addArgument("filter_zero", "", mitkCommandLineParser::Bool, "Filter Zero Weights:", "filter fibers with zero weight", false); parser.addArgument("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("join_tracts", "", mitkCommandLineParser::Bool, "Join output tracts:", "outout tracts are merged into a single tractogram", false); parser.addArgument("regu", "", mitkCommandLineParser::String, "Regularization:", "MSM; Variance; VoxelVariance; Lasso; GroupLasso; GroupVariance; NONE", std::string("VoxelVariance")); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType fib_files = us::any_cast(parsedArgs["i1"]); std::string input_image_name = us::any_cast(parsedArgs["i2"]); std::string outRoot = us::any_cast(parsedArgs["o"]); bool single_fib = true; if (parsedArgs.count("bundle_based")) single_fib = !us::any_cast(parsedArgs["bundle_based"]); bool save_residuals = false; if (parsedArgs.count("save_res")) save_residuals = us::any_cast(parsedArgs["save_res"]); bool filter_zero = false; if (parsedArgs.count("filter_zero")) filter_zero = us::any_cast(parsedArgs["filter_zero"]); bool save_weights = false; if (parsedArgs.count("save_weights")) save_weights = us::any_cast(parsedArgs["save_weights"]); std::string regu = "VoxelVariance"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool join_tracts = false; if (parsedArgs.count("join_tracts")) join_tracts = us::any_cast(parsedArgs["join_tracts"]); int max_iter = 20; if (parsedArgs.count("max_iter")) max_iter = us::any_cast(parsedArgs["max_iter"]); float g_tol = 1e-5; if (parsedArgs.count("min_g")) g_tol = us::any_cast(parsedArgs["min_g"]); float lambda = 1.0; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); try { MITK_INFO << "Loading data"; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {}); std::vector< std::string > fib_names; auto input_tracts = mitk::DiffusionDataIOHelper::load_fibs(fib_files, &fib_names); itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input_image_name, &functor)[0].GetPointer(); mitk::Image::Pointer mitk_image = dynamic_cast(inputData.GetPointer()); mitk::PeakImage::Pointer mitk_peak_image = dynamic_cast(inputData.GetPointer()); if (mitk_peak_image.IsNotNull()) { typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitk_peak_image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer peak_image = caster->GetOutput(); fitter->SetPeakImage(peak_image); } else if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { fitter->SetDiffImage(mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_image)); mitk::TensorModel<>* model = new mitk::TensorModel<>(); model->SetBvalue(1000); model->SetDiffusivity1(0.0010); model->SetDiffusivity2(0.00015); model->SetDiffusivity3(0.00015); model->SetGradientList(mitk::DiffusionPropertyHelper::GetGradientContainer(mitk_image)); fitter->SetSignalModel(model); } else if (mitk_image->GetDimension()==3) { itk::FitFibersToImageFilter::DoubleImgType::Pointer scalar_image = itk::FitFibersToImageFilter::DoubleImgType::New(); mitk::CastToItkImage(mitk_image, scalar_image); fitter->SetScalarImage(scalar_image); } else { MITK_INFO << "Input image invalid. Valid options are peak image, 3D scalar image or raw diffusion-weighted image."; return EXIT_FAILURE; } fitter->SetTractograms(input_tracts); fitter->SetFitIndividualFibers(single_fib); fitter->SetMaxIterations(max_iter); fitter->SetGradientTolerance(g_tol); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); mitk::LocaleSwitch localeSwitch("C"); if (save_residuals && mitk_peak_image.IsNotNull()) { itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New(); writer->SetInput(fitter->GetFittedImage()); writer->SetFileName(outRoot + "_fitted.nii.gz"); writer->Update(); writer->SetInput(fitter->GetResidualImage()); writer->SetFileName(outRoot + "_residual.nii.gz"); writer->Update(); writer->SetInput(fitter->GetOverexplainedImage()); writer->SetFileName(outRoot + "_overexplained.nii.gz"); writer->Update(); writer->SetInput(fitter->GetUnderexplainedImage()); writer->SetFileName(outRoot + "_underexplained.nii.gz"); writer->Update(); } else if (save_residuals && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetFittedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_fitted_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetResidualImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_residual_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetOverexplainedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_overexplained_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetUnderexplainedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_underexplained_image.nii.gz"); } } else if (save_residuals) { itk::ImageFileWriter< itk::FitFibersToImageFilter::DoubleImgType >::Pointer writer = itk::ImageFileWriter< itk::FitFibersToImageFilter::DoubleImgType >::New(); writer->SetInput(fitter->GetFittedImageScalar()); writer->SetFileName(outRoot + "_fitted_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetResidualImageScalar()); writer->SetFileName(outRoot + "_residual_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetOverexplainedImageScalar()); writer->SetFileName(outRoot + "_overexplained_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetUnderexplainedImageScalar()); writer->SetFileName(outRoot + "_underexplained_image.nii.gz"); writer->Update(); } std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); if (!join_tracts) { for (unsigned int bundle=0; bundleFilterByWeights(0.0); if (fib->GetNumFibers()>0) { fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(fib, outRoot + name + "_fitted.fib"); if (save_weights) { ofstream logfile; logfile.open (outRoot + name + "_weights.txt"); for (unsigned int f=0; fGetNumFibers(); ++f) logfile << output_tracts.at(bundle)->GetFiberWeight(f) << "\n"; logfile.close(); } } else MITK_INFO << "Output contains no fibers!"; } } else { mitk::FiberBundle::Pointer out = mitk::FiberBundle::New(); out = out->AddBundles(output_tracts); if (filter_zero) out = out->FilterByWeights(0.0); if (out->GetNumFibers()>0) { out->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out, outRoot + "_fitted.fib"); if (save_weights) { ofstream logfile; logfile.open (outRoot + "_weights.txt"); for (unsigned int f=0; fGetNumFibers(); ++f) logfile << out->GetFiberWeight(f) << "\n"; logfile.close(); } } else MITK_INFO << "Output contains no fibers!"; } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/Sift2WeightCopy.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/Sift2WeightCopy.cpp index c74e325737..e0cd36395c 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/Sift2WeightCopy.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/Sift2WeightCopy.cpp @@ -1,107 +1,107 @@ /*=================================================================== 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::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Import Sift2 Fiber Weights txt file. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("SIFT2 Fiber Weight Import"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Import SIFT2 fiber weights."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output fiber bundle (.fib)", us::Any(), false); parser.addArgument("weights", "", mitkCommandLineParser::String, "Weights:", "input weights file (.txt)", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string weightsFileName = us::any_cast(parsedArgs["weights"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { mitk::FiberBundle::Pointer fib = LoadFib(inFileName); std::ifstream fin; fin.open(weightsFileName); if (!fin.good()) return 1; // exit if file not found std::vector weights; for (float d; fin >> d; ) { weights.push_back(d); } for(std::size_t i = 0; i != weights.size(); i++) { fib->SetFiberWeight(i, weights[i]); } mitk::IOUtil::Save(fib.GetPointer(), outFileName ); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/TractDensity.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensity.cpp index 68e9791c05..0fdc8162da 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensity.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensity.cpp @@ -1,213 +1,213 @@ /*=================================================================== 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::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Modify input tractogram: fiber resampling, compression, pruning and transformation. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Tract Density"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Generate tract density image, fiber envelope or fiber endpoints image."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false); parser.addArgument("binary", "", mitkCommandLineParser::Bool, "Binary output:", "calculate binary tract envelope", us::Any()); parser.addArgument("normalize", "", mitkCommandLineParser::Bool, "Normalized output:", "normalize output to 0-1", us::Any()); parser.addArgument("endpoints", "", mitkCommandLineParser::Bool, "Output endpoints image:", "calculate image of fiber endpoints instead of mask", us::Any()); parser.addArgument("reference_image", "", mitkCommandLineParser::String, "Reference image:", "output image will have geometry of this reference image", us::Any()); parser.addArgument("upsampling", "", mitkCommandLineParser::Float, "Upsampling:", "upsampling", 1.0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; bool binary = false; if (parsedArgs.count("binary")) binary = us::any_cast(parsedArgs["binary"]); bool endpoints = false; if (parsedArgs.count("endpoints")) endpoints = us::any_cast(parsedArgs["endpoints"]); bool normalize = false; if (parsedArgs.count("normalize")) normalize = us::any_cast(parsedArgs["normalize"]); float upsampling = 1.0; if (parsedArgs.count("upsampling")) upsampling = us::any_cast(parsedArgs["upsampling"]); MITK_INFO << "Upsampling: " << upsampling; std::string reference_image = ""; if (parsedArgs.count("reference_image")) reference_image = us::any_cast(parsedArgs["reference_image"]); std::string inFileName = us::any_cast(parsedArgs["i"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { mitk::FiberBundle::Pointer fib = LoadFib(inFileName); mitk::Image::Pointer ref_img; if (!reference_image.empty()) ref_img = mitk::IOUtil::Load(reference_image); if (endpoints) { typedef unsigned int OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(upsampling); if (ref_img.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(ref_img, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, outFileName ); } else if (binary) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(!normalize); generator->SetUpsamplingFactor(upsampling); if (ref_img.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(ref_img, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, outFileName ); } else { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(!normalize); generator->SetUpsamplingFactor(upsampling); if (ref_img.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(ref_img, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, outFileName ); } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp index 4098cc2de9..125cf41a18 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp @@ -1,110 +1,110 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Filter Outliers by Tract Density"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", 0.05); parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap:", "positive means ROI image value threshold", 0.5); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. num. fibers:", "discard positive tracts with less fibers", 0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["i"]); std::string outFib = us::any_cast(parsedArgs["o"]); int min_fibers = 0; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); float overlap = 0.5; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); float threshold = 0.05f; if (parsedArgs.count("threshold")) threshold = us::any_cast(parsedArgs["threshold"]); try { mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(inFib); itk::TractDensityImageFilter< ItkFloatImgType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImgType >::New(); generator->SetFiberBundle(inputTractogram); generator->SetBinaryOutput(false); generator->SetOutputAbsoluteValues(false); generator->Update(); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetRoiImages({generator->GetOutput()}); extractor->SetInputFiberBundle(inputTractogram); extractor->SetOverlapFraction(overlap); extractor->SetInterpolate(true); extractor->SetThreshold(threshold); extractor->SetNoNegatives(true); extractor->Update(); if (extractor->GetPositives().at(0)->GetNumFibers() >= static_cast(min_fibers)) mitk::IOUtil::Save(extractor->GetPositives().at(0), outFib); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Fiberfox/RandomFiberPhantom.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/RandomFiberPhantom.cpp index 11e486d5c4..12d4d29310 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/RandomFiberPhantom.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/RandomFiberPhantom.cpp @@ -1,172 +1,172 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Random Fiber Phantom"); parser.setCategory("Diffusion Simulation Tools"); parser.setContributor("MIC"); parser.setDescription("Create Random Fiber Configurations"); parser.setArgumentPrefix("--", "-"); parser.addArgument("num_bundles", "", mitkCommandLineParser::Int, "", "", 50); parser.addArgument("min_density", "", mitkCommandLineParser::Int, "", "", 50); parser.addArgument("max_density", "", mitkCommandLineParser::Int, "", "", 200); parser.addArgument("size_x", "", mitkCommandLineParser::Int, "", "", 250); parser.addArgument("size_y", "", mitkCommandLineParser::Int, "", "", 250); parser.addArgument("size_z", "", mitkCommandLineParser::Int, "", "", 250); parser.addArgument("min_stepsize", "", mitkCommandLineParser::Int, "", "", 15); parser.addArgument("max_stepsize", "", mitkCommandLineParser::Int, "", "", 30); parser.addArgument("min_curve", "", mitkCommandLineParser::Int, "", "", 5); parser.addArgument("max_curve", "", mitkCommandLineParser::Int, "", "", 45); parser.addArgument("min_radius", "", mitkCommandLineParser::Int, "", "", 5); parser.addArgument("max_radius", "", mitkCommandLineParser::Int, "", "", 25); parser.addArgument("min_twist", "", mitkCommandLineParser::Int, "", "", 15); parser.addArgument("max_twist", "", mitkCommandLineParser::Int, "", "", 30); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "compress fiber using the given error threshold (in mm)", 0.1); parser.addArgument("", "o", mitkCommandLineParser::String, "Output folder:", "output folder", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("fix_seed", "", mitkCommandLineParser::Int, "Fix random seed:", "if >= 0, produce same random values on each run using this seed.", -1); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string out_folder = us::any_cast(parsedArgs["o"]); float compress=0.1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); int num_bundles=50; if (parsedArgs.count("num_bundles")) num_bundles = us::any_cast(parsedArgs["num_bundles"]); int min_density=50; if (parsedArgs.count("min_density")) min_density = us::any_cast(parsedArgs["min_density"]); int max_density=200; if (parsedArgs.count("max_density")) max_density = us::any_cast(parsedArgs["max_density"]); int size_x=250; if (parsedArgs.count("size_x")) size_x = us::any_cast(parsedArgs["size_x"]); int size_y=250; if (parsedArgs.count("size_y")) size_y = us::any_cast(parsedArgs["size_y"]); int size_z=250; if (parsedArgs.count("size_z")) size_z = us::any_cast(parsedArgs["size_z"]); int min_stepsize=15; if (parsedArgs.count("min_stepsize")) min_stepsize = us::any_cast(parsedArgs["min_stepsize"]); int max_stepsize=30; if (parsedArgs.count("max_stepsize")) max_stepsize = us::any_cast(parsedArgs["max_stepsize"]); int min_curve=5; if (parsedArgs.count("min_curve")) min_curve = us::any_cast(parsedArgs["min_curve"]); int max_curve=45; if (parsedArgs.count("max_curve")) max_curve = us::any_cast(parsedArgs["max_curve"]); int min_radius=5; if (parsedArgs.count("min_radius")) min_radius = us::any_cast(parsedArgs["min_radius"]); int max_radius=25; if (parsedArgs.count("max_radius")) max_radius = us::any_cast(parsedArgs["max_radius"]); int min_twist=15; if (parsedArgs.count("min_twist")) min_twist = us::any_cast(parsedArgs["min_twist"]); int max_twist=30; if (parsedArgs.count("max_twist")) max_twist = us::any_cast(parsedArgs["max_twist"]); int fix_seed = -1; if (parsedArgs.count("fix_seed")) fix_seed = us::any_cast(parsedArgs["fix_seed"]); try { itk::RandomPhantomFilter::Pointer filter = itk::RandomPhantomFilter::New(); filter->SetNumTracts(static_cast(num_bundles)); filter->SetMinStreamlineDensity(static_cast(min_density)); filter->SetMaxStreamlineDensity(static_cast(max_density)); mitk::Vector3D vol; vol[0] = size_x; vol[1] = size_y; vol[2] = size_z; filter->SetVolumeSize(vol); filter->SetStepSizeMin(static_cast(min_stepsize)); filter->SetStepSizeMax(static_cast(max_stepsize)); filter->SetCurvynessMin(static_cast(min_curve)); filter->SetCurvynessMax(static_cast(max_curve)); filter->SetStartRadiusMin(static_cast(min_radius)); filter->SetStartRadiusMax(static_cast(max_radius)); filter->SetMinTwist(static_cast(min_twist)); filter->SetMaxTwist(static_cast(max_twist)); filter->SetFixSeed(fix_seed); filter->Update(); auto fibs = filter->GetFiberBundles(); std::vector< mitk::DataNode::Pointer > fiber_nodes; int c = 1; for (auto fib : fibs) { if (compress>0) fib->Compress(compress); mitk::IOUtil::Save(fib, out_folder + "Bundle_" + boost::lexical_cast(c) + ".fib"); ++c; } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/CopyGeometry.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/CopyGeometry.cpp index 1ad8350099..6ea85d3c7d 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/CopyGeometry.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/CopyGeometry.cpp @@ -1,114 +1,114 @@ /*=================================================================== 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 "mitkCommandLineParser.h" /*! \brief Copies transformation matrix of one image to another */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Copy Geometry"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Copies transformation matrix of one image to another"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("ref", "r", mitkCommandLineParser::String, "Reference:", "reference image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("alignCentroid", "a", mitkCommandLineParser::Bool, "align centroids", "align centroids", us::Any(), true); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string refImage = us::any_cast(parsedArgs["ref"]); std::string outImage = us::any_cast(parsedArgs["o"]); bool originOnly = false; // Show a help message if ( parsedArgs.count("alignCentroid") || parsedArgs.count("a")) { originOnly = true; } try { mitk::Image::Pointer source = mitk::IOUtil::Load(refImage); mitk::Image::Pointer target = mitk::IOUtil::Load(imageName); if (originOnly) { // Calculate correction to align centroids double c[3]; c[0] = source->GetGeometry()->GetOrigin()[0] + source->GetGeometry()->GetExtent(0)/2.0 - target->GetGeometry()->GetOrigin()[0] - target->GetGeometry()->GetExtent(0)/2.0; c[1] = source->GetGeometry()->GetOrigin()[1] + source->GetGeometry()->GetExtent(1)/2.0 - target->GetGeometry()->GetOrigin()[1] - target->GetGeometry()->GetExtent(1)/2.0; c[2] = source->GetGeometry()->GetOrigin()[2] + source->GetGeometry()->GetExtent(2)/2.0 - target->GetGeometry()->GetOrigin()[2] - target->GetGeometry()->GetExtent(2)/2.0; double newOrigin[3]; newOrigin[0] = target->GetGeometry()->GetOrigin()[0] +c[0]; newOrigin[1] = target->GetGeometry()->GetOrigin()[1] +c[1]; newOrigin[2] = target->GetGeometry()->GetOrigin()[2] +c[2]; target->GetGeometry()->SetOrigin(newOrigin); } else { mitk::BaseGeometry* s_geom = source->GetGeometry(); mitk::BaseGeometry* t_geom = target->GetGeometry(); t_geom->SetIndexToWorldTransform(s_geom->GetIndexToWorldTransform()); target->SetGeometry(t_geom); } mitk::IOUtil::Save(target, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/DImp.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DImp.cpp index 8b9b9a7a4f..a8b6a06238 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DImp.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DImp.cpp @@ -1,72 +1,72 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("DIMP"); parser.setCategory("Preprocessing Tools"); parser.setDescription("TEMPORARY: Converts DICOM to other image types"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); mitk::Image::Pointer source = mitk::IOUtil::Load(imageName, &functor); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(source, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(source, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/DReg.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DReg.cpp index 1817dfcc95..96dc2ccabe 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DReg.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DReg.cpp @@ -1,222 +1,222 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include typedef mitk::DiffusionPropertyHelper DPH; mitk::Image::Pointer apply_transform(mitk::Image::Pointer moving, mitk::Image::Pointer fixed_single, mitk::MAPRegistrationWrapper::Pointer reg, bool resample) { mitk::Image::Pointer registered_image; if (!resample) { registered_image = mitk::ImageMappingHelper::refineGeometry(moving, reg, true); } else { if (!mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(moving)) { registered_image = mitk::ImageMappingHelper::map(moving, reg, false, 0, fixed_single->GetGeometry(), false, 0, mitk::ImageMappingInterpolator::BSpline_3); } else { typedef itk::Image ITKDiffusionVolumeType; typedef itk::ComposeImageFilter < ITKDiffusionVolumeType > ComposeFilterType; auto composer = ComposeFilterType::New(); auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(moving); for (unsigned int i=0; iGetVectorLength(); ++i) { itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); filter->SetInput( itkVectorImagePointer); filter->SetChannelIndex(i); filter->Update(); mitk::Image::Pointer gradientVolume = mitk::Image::New(); gradientVolume->InitializeByItk( filter->GetOutput() ); gradientVolume->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); mitk::Image::Pointer registered_mitk_image = mitk::ImageMappingHelper::map(gradientVolume, reg, false, 0, fixed_single->GetGeometry(), false, 0, mitk::ImageMappingInterpolator::BSpline_3); auto registered_itk_image = ITKDiffusionVolumeType::New(); mitk::CastToItkImage(registered_mitk_image, registered_itk_image); composer->SetInput(i, registered_itk_image); } composer->Update(); registered_image = mitk::GrabItkImageMemory( composer->GetOutput() ); mitk::DiffusionPropertyHelper::CopyProperties(moving, registered_image, true); typedef mitk::DiffusionImageCorrectionFilter CorrectionFilterType; CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); corrector->SetImage( registered_image ); corrector->CorrectDirections( mitk::MITKRegistrationHelper::getAffineMatrix(reg, false)->GetMatrix().GetVnlMatrix() ); } } return registered_image; } int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("DREG"); parser.setCategory("Preprocessing Tools"); parser.setDescription("TEMPORARY: Rigid registration of two images"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "f", mitkCommandLineParser::String, "Fixed:", "fixed image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "m", mitkCommandLineParser::String, "Moving:", "moving image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("resample", "", mitkCommandLineParser::Bool, "Resample:", "resample moving image", false); parser.addArgument("coreg", "", mitkCommandLineParser::StringList, "", "additionally apply transform to these images", us::Any(), true, false, false, mitkCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string f = us::any_cast(parsedArgs["f"]); std::string m = us::any_cast(parsedArgs["m"]); std::string o = us::any_cast(parsedArgs["o"]); bool resample = false; if (parsedArgs.count("resample")) resample = true; mitkCommandLineParser::StringContainerType coreg; if (parsedArgs.count("coreg")) coreg = us::any_cast(parsedArgs["coreg"]); try { typedef itk::Image< float, 3 > ItkFloatImageType; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); mitk::Image::Pointer fixed = mitk::IOUtil::Load(f, &functor); mitk::Image::Pointer moving = mitk::IOUtil::Load(m, &functor); mitk::Image::Pointer fixed_single = fixed; mitk::Image::Pointer moving_single = moving; mitk::MultiModalRigidDefaultRegistrationAlgorithm< ItkFloatImageType >::Pointer algo = mitk::MultiModalRigidDefaultRegistrationAlgorithm< ItkFloatImageType >::New(); mitk::MITKAlgorithmHelper helper(algo); if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(fixed)) { DPH::ImageType::Pointer itkVectorImagePointer = DPH::ImageType::New(); mitk::CastToItkImage(fixed, itkVectorImagePointer); itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); filter->SetInput( itkVectorImagePointer); filter->SetChannelIndex(0); filter->Update(); fixed_single = mitk::Image::New(); fixed_single->InitializeByItk( filter->GetOutput() ); fixed_single->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); } if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(moving)) { DPH::ImageType::Pointer itkVectorImagePointer = DPH::ImageType::New(); mitk::CastToItkImage(moving, itkVectorImagePointer); itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); filter->SetInput( itkVectorImagePointer); filter->SetChannelIndex(0); filter->Update(); moving_single = mitk::Image::New(); moving_single->InitializeByItk( filter->GetOutput() ); moving_single->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); } helper.SetData(moving_single, fixed_single); mitk::MAPRegistrationWrapper::Pointer reg = helper.GetMITKRegistrationWrapper(); mitk::Image::Pointer registered_image = apply_transform(moving, fixed_single, reg, resample); if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(registered_image)) { mitk::DiffusionPropertyHelper::InitializeImage( registered_image ); std::string file_extension = itksys::SystemTools::GetFilenameExtension(o); if (file_extension==".nii" || file_extension==".nii.gz") mitk::IOUtil::Save(registered_image, "DWI_NIFTI", o); else mitk::IOUtil::Save(registered_image, o); } else mitk::IOUtil::Save(registered_image, o); std::string path = ist::GetFilenamePath(o) + "/"; std::vector< std::string > file_names; auto coreg_images = mitk::DiffusionDataIOHelper::load_mitk_images(coreg, &file_names); for (unsigned int i=0; i #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("DmriDenoising"); parser.setCategory("Preprocessing Tools"); parser.setDescription("dMRI denoising tool"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("type", "", mitkCommandLineParser::Int, "Type:", "0 (TotalVariation), 1 (Gauss), 2 (NLM)", 0); parser.endGroup(); parser.beginGroup("2. Total variation parameters:"); parser.addArgument("tv_iterations", "", mitkCommandLineParser::Int, "Iterations:", "", 1); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "", 0.1); parser.endGroup(); parser.beginGroup("3. Gauss parameters:"); parser.addArgument("variance", "", mitkCommandLineParser::Float, "Variance:", "", 1.0); parser.endGroup(); parser.beginGroup("4. NLM parameters:"); parser.addArgument("nlm_iterations", "", mitkCommandLineParser::Int, "Iterations:", "", 4); parser.addArgument("sampling_radius", "", mitkCommandLineParser::Int, "Sampling radius:", "", 4); parser.addArgument("patch_radius", "", mitkCommandLineParser::Int, "Patch radius:", "", 1); parser.addArgument("num_patches", "", mitkCommandLineParser::Int, "Num. patches:", "", 10); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); int type = 0; if (parsedArgs.count("type")) type = us::any_cast(parsedArgs["type"]); int tv_iterations = 1; if (parsedArgs.count("tv_iterations")) tv_iterations = us::any_cast(parsedArgs["tv_iterations"]); float lambda = 0.1; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); float variance = 1.0; if (parsedArgs.count("variance")) variance = us::any_cast(parsedArgs["variance"]); int nlm_iterations = 4; if (parsedArgs.count("nlm_iterations")) nlm_iterations = us::any_cast(parsedArgs["nlm_iterations"]); int sampling_radius = 4; if (parsedArgs.count("sampling_radius")) sampling_radius = us::any_cast(parsedArgs["sampling_radius"]); int patch_radius = 1; if (parsedArgs.count("patch_radius")) patch_radius = us::any_cast(parsedArgs["patch_radius"]); int num_patches = 10; if (parsedArgs.count("num_patches")) num_patches = us::any_cast(parsedArgs["num_patches"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); mitk::Image::Pointer input_image = mitk::IOUtil::Load(imageName, &functor); typedef short DiffusionPixelType; typedef itk::VectorImage DwiImageType; typedef itk::Image DwiVolumeType; typedef itk::DiscreteGaussianImageFilter < DwiVolumeType, DwiVolumeType > GaussianFilterType; typedef itk::PatchBasedDenoisingImageFilter < DwiVolumeType, DwiVolumeType > NlmFilterType; typedef itk::VectorImageToImageFilter < DiffusionPixelType > ExtractFilterType; typedef itk::ComposeImageFilter < itk::Image > ComposeFilterType; if (!mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input_image)) mitkThrow() << "Input is not a diffusion-weighted image!"; DwiImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input_image); mitk::Image::Pointer denoised_image = nullptr; switch (type) { case 0: { ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i=0; iGetVectorLength(); ++i) { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); extractor->SetIndex( i ); extractor->Update(); DwiVolumeType::Pointer gradient_volume = extractor->GetOutput(); itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::Pointer filter = itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::New(); filter->SetInput(gradient_volume); filter->SetLambda(lambda); filter->SetNumberIterations(tv_iterations); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); break; } case 1: { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); GaussianFilterType::Pointer filter = GaussianFilterType::New(); filter->SetVariance(variance); filter->SetInput(extractor->GetOutput()); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); break; } case 2: { typedef itk::Statistics::GaussianRandomSpatialNeighborSubsampler< NlmFilterType::PatchSampleType, DwiVolumeType::RegionType > SamplerType; // sampling the image to find similar patches SamplerType::Pointer sampler = SamplerType::New(); sampler->SetRadius( sampling_radius ); sampler->SetVariance( sampling_radius*sampling_radius ); sampler->SetNumberOfResultsRequested( num_patches ); MITK_INFO << "Starting NLM denoising"; ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); NlmFilterType::Pointer filter = NlmFilterType::New(); filter->SetInput(extractor->GetOutput()); filter->SetPatchRadius(patch_radius); filter->SetNoiseModel(NlmFilterType::RICIAN); filter->UseSmoothDiscPatchWeightsOn(); filter->UseFastTensorComputationsOn(); filter->SetNumberOfIterations(nlm_iterations); filter->SetSmoothingWeight( 1 ); filter->SetKernelBandwidthEstimation(true); filter->SetSampler( sampler ); filter->Update(); composer->SetInput(i, filter->GetOutput()); MITK_INFO << "Gradient " << i << " finished"; } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); break; } } mitk::DiffusionPropertyHelper::SetGradientContainer(denoised_image, mitk::DiffusionPropertyHelper::GetGradientContainer(input_image)); mitk::DiffusionPropertyHelper::SetReferenceBValue(denoised_image, mitk::DiffusionPropertyHelper::GetReferenceBValue(input_image)); mitk::DiffusionPropertyHelper::InitializeImage( denoised_image ); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(denoised_image, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(denoised_image, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/FileFormatConverter.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FileFormatConverter.cpp index a4e0633ba1..2e6eea39b0 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FileFormatConverter.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FileFormatConverter.cpp @@ -1,80 +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 "mitkCommandLineParser.h" #include #include /*! \brief Load image and save as specified file type. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Format Converter"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Load image and save as specified file type."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output file", us::Any(), false, false, false, mitkCommandLineParser::Output); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string inName = us::any_cast(parsedArgs["i"]); std::string outName = us::any_cast(parsedArgs["o"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); std::vector baseData = mitk::IOUtil::Load(inName, &functor); if ( baseData.size()>0 && dynamic_cast(baseData[0].GetPointer()) ) { mitk::IOUtil::Save(dynamic_cast(baseData[0].GetPointer()), outName.c_str()); } else if ( baseData.size()>0 && dynamic_cast(baseData[0].GetPointer()) ) { mitk::IOUtil::Save(dynamic_cast(baseData[0].GetPointer()) ,outName.c_str()); } else std::cout << "File type currently not supported!"; } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/FlipPeaks.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FlipPeaks.cpp index eb9ebaee8c..758265b0e9 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FlipPeaks.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FlipPeaks.cpp @@ -1,104 +1,104 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include /*! \brief Copies transformation matrix of one image to another */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Flip Peaks"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Flips the peaks of the input peak image along the given dimensions."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("", "x", mitkCommandLineParser::Bool, "Flip x", "flip along x-axis"); parser.addArgument("", "y", mitkCommandLineParser::Bool, "Flip y", "flip along y-axis"); parser.addArgument("", "z", mitkCommandLineParser::Bool, "Flip z", "flip along z-axis"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); bool x = false; if (parsedArgs.count("x")) x = us::any_cast(parsedArgs["x"]); bool y = false; if (parsedArgs.count("y")) y = us::any_cast(parsedArgs["y"]); bool z = false; if (parsedArgs.count("z")) z = us::any_cast(parsedArgs["z"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); mitk::PeakImage::Pointer image = mitk::IOUtil::Load(imageName, &functor); typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer itkImg = caster->GetOutput(); itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); flipper->SetInput(itkImg); flipper->SetFlipX(x); flipper->SetFlipY(y); flipper->SetFlipZ(z); flipper->Update(); mitk::Image::Pointer resultImage = dynamic_cast(mitk::PeakImage::New().GetPointer()); mitk::CastToMitkImage(flipper->GetOutput(), resultImage); resultImage->SetVolume(flipper->GetOutput()->GetBufferPointer()); mitk::IOUtil::Save(resultImage, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp index b434a943bf..6633946b93 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp @@ -1,78 +1,78 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("HeadMotionCorrection"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Simple affine head-motion correction tool"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); mitk::Image::Pointer in_image = mitk::IOUtil::Load(imageName, &functor); mitk::DWIHeadMotionCorrectionFilter::Pointer registerer = mitk::DWIHeadMotionCorrectionFilter::New(); registerer->SetInput(in_image); registerer->Update(); mitk::Image::Pointer out_image = registerer->GetCorrectedImage(); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(out_image, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(out_image, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/PeakExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp index 4134069cc7..5f6d486958 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp @@ -1,294 +1,294 @@ /*=================================================================== 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 #include #include #include template int StartPeakExtraction(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input image", "sh coefficient image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output directory", "output root", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("mask", "", mitkCommandLineParser::String, "Mask", "mask image", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("normalization", "", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); parser.addArgument("numpeaks", "", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); parser.addArgument("rel_peakthr", "", mitkCommandLineParser::Float, "Relative peak threshold", "peak threshold relative to largest peak", 0.4, true); parser.addArgument("abs_peakthr", "", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak magnitude threshold", 0.03, true); parser.addArgument("angular_thr", "", mitkCommandLineParser::Float, "Angular threshold", "in degree", 15); parser.addArgument("shConvention", "", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MRtrix, FSL)", std::string("MRtrix"), true); parser.addArgument("flipX", "", mitkCommandLineParser::Bool, "Flip X", "Flip peaks in x direction"); parser.addArgument("flipY", "", mitkCommandLineParser::Bool, "Flip Y", "Flip peaks in y direction"); parser.addArgument("flipZ", "", mitkCommandLineParser::Bool, "Flip Z", "Flip peaks in z direction"); parser.addArgument("scale_by_gfa", "", mitkCommandLineParser::Bool, "Scale by GFA", "Scale ODF values and peaks by GFA"); parser.setCategory("Preprocessing Tools"); parser.setTitle("Peak Extraction"); parser.setDescription(""); parser.setContributor("MIC"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outRoot = us::any_cast(parsedArgs["o"]); // optional arguments std::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 rel_peakthr = 0.4; if (parsedArgs.count("rel_peakthr")) rel_peakthr = us::any_cast(parsedArgs["rel_peakthr"]); float abs_peakthr = 0.03; if (parsedArgs.count("abs_peakthr")) abs_peakthr = us::any_cast(parsedArgs["abs_peakthr"]); float angular_thr = 15; if (parsedArgs.count("angular_thr")) angular_thr = us::any_cast(parsedArgs["angular_thr"]); angular_thr = cos((float)angular_thr*itk::Math::pi / 180); bool scale_by_gfa = false; if (parsedArgs.count("scale_by_gfa")) scale_by_gfa = us::any_cast(parsedArgs["scale_by_gfa"]); bool flipX = false; if (parsedArgs.count("flipX")) flipX = us::any_cast(parsedArgs["flipX"]); bool flipY = false; if (parsedArgs.count("flipY")) flipY = us::any_cast(parsedArgs["flipY"]); bool flipZ = false; if (parsedArgs.count("flipZ")) flipZ = us::any_cast(parsedArgs["flipZ"]); 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: " << rel_peakthr; std::cout << "abspeakthres: " << abs_peakthr; std::cout << "shOrder: " << shOrder; try { mitk::Image::Pointer image = mitk::IOUtil::Load(imageName); mitk::Image::Pointer mask = mitk::IOUtil::Load(maskImageName); typedef itk::Image ItkUcharImgType; typedef itk::OdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer peak_extraction_filter = MaximaExtractionFilterType::New(); ItkUcharImgType::Pointer itkMaskImage = nullptr; if (mask.IsNotNull()) { try{ itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(mask, itkMaskImage); peak_extraction_filter->SetMaskImage(itkMaskImage); } catch(...) { } } if (parsedArgs.count("shConvention")) { std::string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); if ( convention=="FSL" ) peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::FSL); else peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); } else peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); try{ typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); peak_extraction_filter->SetInput(caster->GetOutput()); } catch(...) { std::cout << "wrong image type"; return EXIT_FAILURE; } peak_extraction_filter->SetMaxNumPeaks(numPeaks); peak_extraction_filter->SetRelativePeakThreshold(rel_peakthr); peak_extraction_filter->SetAbsolutePeakThreshold(abs_peakthr); peak_extraction_filter->SetAngularThreshold(angular_thr); peak_extraction_filter->SetFlipX(flipX); peak_extraction_filter->SetFlipY(flipY); peak_extraction_filter->SetFlipZ(flipZ); peak_extraction_filter->SetScaleByGfa(scale_by_gfa); switch (normalization) { case 0: peak_extraction_filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: peak_extraction_filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: peak_extraction_filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } std::cout << "Starting extraction"; peak_extraction_filter->Update(); mitk::LocaleSwitch localeSwitch("C"); // write direction image { typename MaximaExtractionFilterType::PeakImageType::Pointer itkImg = peak_extraction_filter->GetPeakImage(); std::string outfilename = outRoot; outfilename.append("_PEAKS.nrrd"); typedef itk::ImageFileWriter< typename MaximaExtractionFilterType::PeakImageType > WriterType; typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(itkImg); writer->Update(); } // write num directions image { ItkUcharImgType::Pointer numDirImage = peak_extraction_filter->GetNumDirectionsImage(); if (itkMaskImage.IsNotNull()) { numDirImage->SetDirection(itkMaskImage->GetDirection()); numDirImage->SetOrigin(itkMaskImage->GetOrigin()); } std::string outfilename = outRoot.c_str(); outfilename.append("_NUM_PEAKS.nrrd"); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(numDirImage); writer->Update(); } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } /*! \brief Extract maxima in the input spherical harmonics image. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input image", "sh coefficient image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output directory", "output root", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("shOrder", "sh", mitkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order"); parser.addArgument("mask", "m", mitkCommandLineParser::String, "Mask", "mask image", us::Any(), true, false, false, mitkCommandLineParser::Input); 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)", std::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("Extract maxima in the input spherical harmonics image."); parser.setContributor("MIC"); std::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/DiffusionCmdApps/Misc/ResampleGradients.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ResampleGradients.cpp index 31726f23c3..8d53507b0a 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ResampleGradients.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ResampleGradients.cpp @@ -1,230 +1,230 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include "itkDWIVoxelFunctor.h" #include typedef short DiffusionPixelType; typedef itk::VectorImage< short, 3 > ItkDwiType; // itk includes #include "itkTimeProbe.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "itkBrainMaskExtractionImageFilter.h" #include "itkCastImageFilter.h" #include "itkVectorContainer.h" #include #include #include #include #include #include // Multishell includes #include // Multishell Functors #include #include #include #include // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" //#include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include #include #include #include #include //#include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkPreferenceListReaderOptionsFunctor.h" mitk::Image::Pointer DoReduceGradientDirections(mitk::Image::Pointer image, double BValue, unsigned int numOfGradientsToKeep, bool use_first_n) { bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(image) ); if ( !isDiffusionImage ) { std::cout << "Image is not a Diffusion Weighted Image" << std::endl; //return; } typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; typedef mitk::BValueMapProperty::BValueMap BValueMap; BValueMap shellSlectionMap; BValueMap originalShellMap = mitk::DiffusionPropertyHelper::GetBValueMap(image); std::vector newNumGradientDirections; //Keeps 1 b0 gradient double B0Value = 0; shellSlectionMap[B0Value] = originalShellMap[B0Value]; unsigned int num = 1; newNumGradientDirections.push_back(num); //BValue = 1000; shellSlectionMap[BValue] = originalShellMap[BValue]; //numOfGradientsToKeep = 32; newNumGradientDirections.push_back(numOfGradientsToKeep); if (newNumGradientDirections.empty()) { std::cout << "newNumGradientDirections is empty" << std::endl; //return; } auto gradientContainer = mitk::DiffusionPropertyHelper::GetGradientContainer(image); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(image, itkVectorImagePointer); std::cout << "1" << std::endl; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkVectorImagePointer ); filter->SetOriginalGradientDirections(gradientContainer); filter->SetNumGradientDirections(newNumGradientDirections); filter->SetOriginalBValueMap(originalShellMap); filter->SetShellSelectionBValueMap(shellSlectionMap); filter->SetUseFirstN(use_first_n); filter->Update(); std::cout << "2" << std::endl; if( filter->GetOutput() == nullptr) { std::cout << "filter get output is nullptr" << std::endl; } mitk::Image::Pointer newImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::CopyProperties(image, newImage, true); mitk::DiffusionPropertyHelper::SetGradientContainer(newImage, filter->GetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( newImage ); return newImage; } /*! \brief Resample gradients of input DWI image. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Resample Gradients"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Resample gradients of input DWI image. You can select one b-value shell and the number of gradients within this shell you want to have. It will also keep one b0 image."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false); parser.addArgument("b_value", "", mitkCommandLineParser::Float, "b-value:", "float", 1000, false); parser.addArgument("num_gradients", "", mitkCommandLineParser::Int, "Nr of gradients:", "integer", 32, false); parser.addArgument("use_first_n", "", mitkCommandLineParser::Bool, "Use first N:", "no optimization, simply use first n gradients", 0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string outFileName = us::any_cast(parsedArgs["o"]); double bValue = us::any_cast(parsedArgs["b_value"]); unsigned int nrOfGradients = us::any_cast(parsedArgs["num_gradients"]); bool use_first_n = false; if (parsedArgs.count("use_first_n")) use_first_n = true; try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({ "Diffusion Weighted Images" }, {}); mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(inFileName, &functor); mitk::Image::Pointer newImage = DoReduceGradientDirections(mitkImage, bValue, nrOfGradients, use_first_n); //mitk::IOUtil::Save(newImage, outFileName); //save as dwi image mitk::IOUtil::Save(newImage, "DWI_NIFTI", outFileName); //save as nifti image } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/RoundBvalues.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/RoundBvalues.cpp index 966a1384d6..c7170b97e8 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/RoundBvalues.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/RoundBvalues.cpp @@ -1,106 +1,106 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("RoundBvalues"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Round b-values"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("to_nearest", "", mitkCommandLineParser::Int, "To nearest:", "integer", 1000); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); int to_nearest = 1000; if (parsedArgs.count("to_nearest")) to_nearest = us::any_cast(parsedArgs["to_nearest"]); try { typedef mitk::DiffusionPropertyHelper PropHelper; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); mitk::Image::Pointer in_image = mitk::IOUtil::Load(imageName, &functor); if (!PropHelper::IsDiffusionWeightedImage(in_image)) { mitkThrow() << "Input is not a diffusion weighted image: " << imageName; } typedef itk::DwiGradientLengthCorrectionFilter FilterType; auto itkVectorImagePointer = PropHelper::GetItkVectorImage(in_image); FilterType::Pointer filter = FilterType::New(); filter->SetRoundingValue(to_nearest); filter->SetReferenceBValue(PropHelper::GetReferenceBValue(in_image)); filter->SetReferenceGradientDirectionContainer(PropHelper::GetGradientContainer(in_image)); filter->Update(); mitk::Image::Pointer newImage = mitk::Image::New(); newImage->InitializeByItk( itkVectorImagePointer.GetPointer() ); newImage->SetImportVolume( itkVectorImagePointer->GetBufferPointer(), 0, 0, mitk::Image::CopyMemory); itkVectorImagePointer->GetPixelContainer()->ContainerManageMemoryOff(); PropHelper::CopyProperties(in_image, newImage, true); PropHelper::SetReferenceBValue(newImage, filter->GetNewBValue()); PropHelper::SetGradientContainer(newImage, filter->GetOutputGradientDirectionContainer()); PropHelper::InitializeImage(newImage); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(newImage, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(newImage, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Misc/ShToOdfImage.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp index bdd2381160..55b2301be6 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp @@ -1,72 +1,72 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("ShToOdfImage"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Calculate discrete ODF image from SH coefficient image"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); mitk::ShImage::Pointer source = mitk::IOUtil::Load(imageName, &functor); mitk::Image::Pointer mitkImage = dynamic_cast(source.GetPointer()); mitk::OdfImage::Pointer out_image = mitk::convert::GetOdfFromShImage(mitkImage); if (out_image.IsNotNull()) mitk::IOUtil::Save(out_image, outImage); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Quantification/DiffusionIndices.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionIndices.cpp index 510144b159..600874d58e 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionIndices.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionIndices.cpp @@ -1,193 +1,193 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #include /** * */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Diffusion Indices"); parser.setCategory("Diffusion Related Measures"); parser.setDescription("Computes requested diffusion related measures"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("index", "idx", mitkCommandLineParser::String, "Index:", "index (fa, gfa, ra, ad, rd, ca, l2, l3, md, adc)", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string index = us::any_cast(parsedArgs["index"]); std::string outFileName = us::any_cast(parsedArgs["o"]); std::string ext = itksys::SystemTools::GetFilenameLastExtension(outFileName); if (ext.empty()) outFileName += ".nii.gz"; try { // load input image mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "SH Image", "ODF Image", "Tensor Image"}, {}); auto input = mitk::IOUtil::Load(inFileName, &functor); bool is_odf = (dynamic_cast(input.GetPointer()) || dynamic_cast(input.GetPointer())); bool is_dt = dynamic_cast(input.GetPointer()); bool is_dw = mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input); if (is_odf) MITK_INFO << "Input is ODF image"; else if (is_dt) MITK_INFO << "Input is tensor image"; else if (is_dw) MITK_INFO << "Input is dMRI"; else { MITK_WARN << "Input is no ODF, SH, tensor or raw dMRI."; return EXIT_FAILURE; } mitk::LocaleSwitch localeSwitch("C"); if( is_odf && index=="gfa" ) { typedef itk::Vector OdfVectorType; typedef itk::Image OdfVectorImgType; OdfVectorImgType::Pointer itkvol; if (dynamic_cast(input.GetPointer())) itkvol = mitk::convert::GetItkOdfFromShImage(input); else itkvol = mitk::convert::GetItkOdfFromOdfImage(input); typedef itk::DiffusionOdfGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(itkvol); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(gfaFilter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if( is_dt ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(input.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itk_dti.GetPointer() ); if(index=="fa") measurementsCalculator->SetMeasure(MeasurementsType::FA); else if(index=="ra") measurementsCalculator->SetMeasure(MeasurementsType::RA); else if(index=="ad") measurementsCalculator->SetMeasure(MeasurementsType::AD); else if(index=="rd") measurementsCalculator->SetMeasure(MeasurementsType::RD); else if(index=="ca") measurementsCalculator->SetMeasure(MeasurementsType::CA); else if(index=="l2") measurementsCalculator->SetMeasure(MeasurementsType::L2); else if(index=="l3") measurementsCalculator->SetMeasure(MeasurementsType::L3); else if(index=="md") measurementsCalculator->SetMeasure(MeasurementsType::MD); else { MITK_WARN << "No valid diffusion index for input image (tensor image) defined"; return EXIT_FAILURE; } measurementsCalculator->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(measurementsCalculator->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if(is_dw && (index=="adc" || index=="md")) { typedef itk::AdcImageFilter< short, double > FilterType; auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkVectorImagePointer ); filter->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(input) ); filter->SetB_value( static_cast(mitk::DiffusionPropertyHelper::GetReferenceBValue(input)) ); if (index=="adc") filter->SetFitSignal(true); else filter->SetFitSignal(false); filter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(filter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else std::cout << "Diffusion index " << index << " not supported for supplied file type."; } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Quantification/MultishellMethods.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/MultishellMethods.cpp index bab85deda6..63eb7ec37e 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/MultishellMethods.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/MultishellMethods.cpp @@ -1,215 +1,215 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #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("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output file", us::Any(), false, false, false, mitkCommandLineParser::Output); 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); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string inName = us::any_cast(parsedArgs["i"]); std::string outName = us::any_cast(parsedArgs["o"]); bool applyADC = us::any_cast(parsedArgs["adc"]); bool applyAKC = us::any_cast(parsedArgs["akc"]); bool applyBiExp = us::any_cast(parsedArgs["biexp"]); std::string targetType = us::any_cast(parsedArgs["targetbvalue"]); try { std::cout << "Loading " << inName; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); mitk::Image::Pointer dwi = mitk::IOUtil::Load(inName, &functor); if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dwi ) ) { typedef itk::RadialMultishellToSingleshellImageFilter FilterType; typedef itk::DwiGradientLengthCorrectionFilter CorrectionFilterType; CorrectionFilterType::Pointer roundfilter = CorrectionFilterType::New(); roundfilter->SetRoundingValue( 1000 ); roundfilter->SetReferenceBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi )); roundfilter->SetReferenceGradientDirectionContainer(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi)); roundfilter->Update(); mitk::DiffusionPropertyHelper::SetReferenceBValue(dwi, roundfilter->GetNewBValue()); mitk::DiffusionPropertyHelper::SetGradientContainer(dwi, roundfilter->GetOutputGradientDirectionContainer()); // filter input parameter const mitk::DiffusionPropertyHelper::BValueMapType &originalShellMap = mitk::DiffusionPropertyHelper::GetBValueMap(dwi); mitk::DiffusionPropertyHelper::ImageType::Pointer vectorImage = mitk::DiffusionPropertyHelper::ImageType::New(); mitk::CastToItkImage(dwi, vectorImage); const mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainer = mitk::DiffusionPropertyHelper::GetGradientContainer(dwi); const unsigned int &bValue = mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi ); // filter call vnl_vector bValueList(originalShellMap.size()-1); double targetBValue = bValueList.mean(); mitk::DiffusionPropertyHelper::BValueMapType::const_iterator it = originalShellMap.begin(); ++it; int i = 0 ; for(; it != originalShellMap.end(); ++it) bValueList.put(i++,it->first); if( targetType == "mean" ) targetBValue = bValueList.mean(); else if( targetType == "min" ) targetBValue = bValueList.min_value(); else if( targetType == "max" ) targetBValue = bValueList.max_value(); if(applyADC) { FilterType::Pointer filter = FilterType::New(); filter->SetInput(vectorImage); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetOriginalBValue(bValue); itk::ADCAverageFunctor::Pointer functor = itk::ADCAverageFunctor::New(); functor->setListOfBValues(bValueList); functor->setTargetBValue(targetBValue); filter->SetFunctor(functor); filter->Update(); // create new DWI image mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetReferenceBValue(outImage, targetBValue); mitk::DiffusionPropertyHelper::SetGradientContainer(outImage, filter->GetTargetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, (outName + "_ADC.dwi").c_str()); } if(applyAKC) { FilterType::Pointer filter = FilterType::New(); filter->SetInput(vectorImage); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetOriginalBValue(bValue); itk::KurtosisFitFunctor::Pointer functor = itk::KurtosisFitFunctor::New(); functor->setListOfBValues(bValueList); functor->setTargetBValue(targetBValue); filter->SetFunctor(functor); filter->Update(); // create new DWI image mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetReferenceBValue(outImage, targetBValue); mitk::DiffusionPropertyHelper::SetGradientContainer(outImage, filter->GetTargetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, (std::string(outName) + "_AKC.dwi").c_str()); } if(applyBiExp) { FilterType::Pointer filter = FilterType::New(); filter->SetInput(vectorImage); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetOriginalBValue(bValue); itk::BiExpFitFunctor::Pointer functor = itk::BiExpFitFunctor::New(); functor->setListOfBValues(bValueList); functor->setTargetBValue(targetBValue); filter->SetFunctor(functor); filter->Update(); // create new DWI image mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetReferenceBValue(outImage, targetBValue); mitk::DiffusionPropertyHelper::SetGradientContainer(outImage, filter->GetTargetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, (std::string(outName) + "_BiExp.dwi").c_str()); } } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/Tractography/GlobalTractography.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp index 72978ab5d9..61973b7e12 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp @@ -1,133 +1,133 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include /*! \brief Perform global fiber tractography (Gibbs tractography) */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Gibbs Tracking"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform global fiber tractography (Gibbs tractography)"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("parameters", "", mitkCommandLineParser::String, "Parameters:", "parameter file (.gtp)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("mask", "", mitkCommandLineParser::String, "Mask:", "binary mask image", us::Any(), false, false, false, mitkCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string paramFileName = us::any_cast(parsedArgs["parameters"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { // instantiate gibbs tracker typedef itk::Vector OdfVectorType; typedef itk::Image ItkOdfImageType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); // load input image mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(inFileName, &functor); // try to cast to Odf image if( dynamic_cast(mitkImage.GetPointer()) ) { mitk::OdfImage::Pointer mitkOdfImage = dynamic_cast(mitkImage.GetPointer()); ItkOdfImageType::Pointer itk_odf = ItkOdfImageType::New(); mitk::CastToItkImage(mitkOdfImage, itk_odf); gibbsTracker->SetOdfImage(itk_odf.GetPointer()); } else if( dynamic_cast(mitkImage.GetPointer()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(mitkImage.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); gibbsTracker->SetTensorImage(itk_dti); } else if ( dynamic_cast(mitkImage.GetPointer()) ) { mitk::Image::Pointer shImage = dynamic_cast(mitkImage.GetPointer()); gibbsTracker->SetOdfImage(mitk::convert::GetItkOdfFromShImage(shImage)); } else return EXIT_FAILURE; // global tracking if (parsedArgs.count("mask")) { typedef itk::Image MaskImgType; mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::Load(us::any_cast(parsedArgs["mask"])); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); gibbsTracker->SetMaskImage(itk_mask); } gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetLoadParameterFile( paramFileName ); // gibbsTracker->SetLutPath( "" ); gibbsTracker->Update(); mitk::FiberBundle::Pointer mitkFiberBundle = mitk::FiberBundle::New(gibbsTracker->GetFiberBundle()); mitkFiberBundle->SetReferenceGeometry(mitkImage->GetGeometry()); mitk::IOUtil::Save(mitkFiberBundle, outFileName ); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp index cc296bd6db..e8fa990b7e 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp @@ -1,447 +1,447 @@ /*=================================================================== 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 typedef itk::Point PointType4; typedef mitk::PeakImage::ItkPeakImageType PeakImgType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; /*! \brief Score input candidate tracts using ACP analysis */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Anchor Constrained Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription("Score input candidate tracts using ACP analysis"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "a", mitkCommandLineParser::String, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("", "p", mitkCommandLineParser::String, "Input peaks:", "input peak image", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "c", mitkCommandLineParser::StringList, "Candidates:", "Folder(s) or file list of candidate tracts", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output folder:", "output folder", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("reference_mask_folders", "", mitkCommandLineParser::StringList, "Reference Mask Folder(s):", "Folder(s) or file list containing reference tract masks for accuracy evaluation", true, false, false, mitkCommandLineParser::Input); parser.addArgument("reference_peaks_folders", "", mitkCommandLineParser::StringList, "Reference Peaks Folder(s):", "Folder(s) or file list containing reference peak images for accuracy evaluation", true, false, false, mitkCommandLineParser::Input); parser.addArgument("mask", "", mitkCommandLineParser::String, "Mask image:", "scoring is only performed inside the mask image", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("greedy_add", "", mitkCommandLineParser::Bool, "Greedy:", "if enabled, the candidate tracts are not jointly fitted to the residual image but one after the other employing a greedy scheme", false); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); parser.addArgument("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("regu", "", mitkCommandLineParser::String, "Regularization:", "MSM; Variance; VoxelVariance; Lasso; GroupLasso; GroupVariance; NONE", std::string("NONE")); parser.addArgument("use_num_streamlines", "", mitkCommandLineParser::Bool, "Use number of streamlines as score:", "Don't fit candidates, simply use number of streamlines per candidate as score", false); parser.addArgument("use_weights", "", mitkCommandLineParser::Bool, "Use input weights as score:", "Don't fit candidates, simply use first input streamline weight per candidate as score", false); parser.addArgument("filter_zero_weights", "", mitkCommandLineParser::Bool, "Filter zero-weights", "Remove streamlines with weight 0 from candidates", false); parser.addArgument("flipx", "", mitkCommandLineParser::Bool, "Flip x", "flip along x-axis", false); parser.addArgument("flipy", "", mitkCommandLineParser::Bool, "Flip y", "flip along y-axis", false); parser.addArgument("flipz", "", mitkCommandLineParser::Bool, "Flip z", "flip along z-axis", false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string peak_file_name = us::any_cast(parsedArgs["p"]); std::string out_folder = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType candidate_tract_folders = us::any_cast(parsedArgs["c"]); if (!out_folder.empty() && out_folder.back() != '/') out_folder += "/"; bool greedy_add = false; if (parsedArgs.count("greedy_add")) greedy_add = us::any_cast(parsedArgs["greedy_add"]); float lambda = 0.1f; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); bool filter_zero_weights = false; if (parsedArgs.count("filter_zero_weights")) filter_zero_weights = us::any_cast(parsedArgs["filter_zero_weights"]); std::string mask_file = ""; if (parsedArgs.count("mask")) mask_file = us::any_cast(parsedArgs["mask"]); mitkCommandLineParser::StringContainerType reference_mask_files_folders; if (parsedArgs.count("reference_mask_folders")) reference_mask_files_folders = us::any_cast(parsedArgs["reference_mask_folders"]); mitkCommandLineParser::StringContainerType reference_peaks_files_folders; if (parsedArgs.count("reference_peaks_folders")) reference_peaks_files_folders = us::any_cast(parsedArgs["reference_peaks_folders"]); std::string regu = "NONE"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool use_weights = false; if (parsedArgs.count("use_weights")) use_weights = us::any_cast(parsedArgs["use_weights"]); bool use_num_streamlines = false; if (parsedArgs.count("use_num_streamlines")) use_num_streamlines = us::any_cast(parsedArgs["use_num_streamlines"]); bool flipx = false; if (parsedArgs.count("flipx")) flipx = us::any_cast(parsedArgs["flipx"]); bool flipy = false; if (parsedArgs.count("flipy")) flipy = us::any_cast(parsedArgs["flipy"]); bool flipz = false; if (parsedArgs.count("flipz")) flipz = us::any_cast(parsedArgs["flipz"]); try { itk::TimeProbe clock; clock.Start(); if (!ist::PathExists(out_folder)) { MITK_INFO << "Creating output directory"; ist::MakeDirectory(out_folder); } MITK_INFO << "Loading data"; // Load mask file. Fit is only performed inside the mask MITK_INFO << "Loading mask image"; auto mask = mitk::DiffusionDataIOHelper::load_itk_image(mask_file); // Load masks covering the true positives for evaluation purposes MITK_INFO << "Loading reference peaks and masks"; std::vector< std::string > anchor_mask_files; auto reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files_folders, &anchor_mask_files); auto reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peaks_files_folders); // Load peak image MITK_INFO << "Loading peak image"; auto peak_image = mitk::DiffusionDataIOHelper::load_itk_image(peak_file_name); // Load all candidate tracts MITK_INFO << "Loading candidate tracts"; std::vector< std::string > candidate_tract_files; auto input_candidates = mitk::DiffusionDataIOHelper::load_fibs(candidate_tract_folders, &candidate_tract_files); if (flipx || flipy || flipz) { itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); flipper->SetInput(peak_image); flipper->SetFlipX(flipx); flipper->SetFlipY(flipy); flipper->SetFlipZ(flipz); flipper->Update(); peak_image = flipper->GetOutput(); } mitk::LocaleSwitch localeSwitch("C"); itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); ofstream logfile; logfile.open (out_folder + "scores.txt"); double rmse = 0.0; int iteration = 0; std::string name = "NOANCHOR"; if (parsedArgs.count("a")) { // Load reference tractogram consisting of all known tracts std::string anchors_file = us::any_cast(parsedArgs["a"]); mitk::FiberBundle::Pointer anchor_tractogram = mitk::IOUtil::Load(anchors_file); if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) ) { // Fit known tracts to peak image to obtain underexplained image MITK_INFO << "Fit anchor tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetTractograms({anchor_tractogram}); fitter->SetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetPeakImage(peak_image); fitter->SetVerbose(true); fitter->SetMaskImage(mask); fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); rmse = fitter->GetRMSE(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); name = ist::GetFilenameWithoutExtension(anchors_file); mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0); anchor_tracts->SetFiberColors(255,255,255); mitk::IOUtil::Save(anchor_tracts, out_folder + boost::lexical_cast(static_cast(100000*rms_diff[0])) + "_" + name + ".fib"); logfile << name << " " << setprecision(5) << rms_diff[0] << "\n"; peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_" + name + ".nii.gz"); peak_image_writer->Update(); } } if (use_weights || use_num_streamlines) { MITK_INFO << "Using tract weights as scores"; unsigned int c = 0; for (auto fib : input_candidates) { int mod = 1; float score = 0; if (use_weights) { score = fib->GetFiberWeight(0); mod = 100000; } else if (use_num_streamlines) score = fib->GetNumFibers(); fib->ColorFibersByOrientation(); std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast(static_cast(mod*score)) + "_" + bundle_name + ".fib"); unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } float weight_sum = 0; for (unsigned int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << bundle_name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } } else if (!greedy_add) { MITK_INFO << "Fit candidate tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(true); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); fitter->SetTractograms(input_candidates); fitter->SetFitIndividualFibers(true); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); unsigned int c = 0; for (auto fib : input_candidates) { std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect if (filter_zero_weights) fib = fib->FilterByWeights(0); mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib"); unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } float weight_sum = 0; for (unsigned int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << bundle_name << " " << setprecision(5) << rms_diff[c] << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New(); out_fib = out_fib->AddBundles(input_candidates); out_fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib"); peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_AllCandidates.nii.gz"); peak_image_writer->Update(); } else { MITK_INFO << "RMSE: " << setprecision(5) << rmse; // fitter->SetPeakImage(peak_image); // Iteratively add candidate bundles in a greedy manner while (!input_candidates.empty()) { double next_rmse = rmse; mitk::FiberBundle::Pointer best_candidate = nullptr; PeakImgType::Pointer best_candidate_peak_image = nullptr; for (unsigned int i=0; iSetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(false); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); // ****************************** fitter->SetTractograms({input_candidates.at(i)}); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect fitter->Update(); std::cout.rdbuf (old); // <-- restore double candidate_rmse = fitter->GetRMSE(); if (candidate_rmseGetTractograms().at(0); best_candidate_peak_image = fitter->GetUnderexplainedImage(); } } if (best_candidate.IsNull()) break; // fitter->SetPeakImage(peak_image); peak_image = best_candidate_peak_image; unsigned int i=0; std::vector< mitk::FiberBundle::Pointer > remaining_candidates; std::vector< std::string > remaining_candidate_files; for (auto fib : input_candidates) { if (fib!=best_candidate) { remaining_candidates.push_back(fib); remaining_candidate_files.push_back(candidate_tract_files.at(i)); } else name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i)); ++i; } input_candidates = remaining_candidates; candidate_tract_files = remaining_candidate_files; iteration++; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect // Save winning candidate if (filter_zero_weights) best_candidate = best_candidate->FilterByWeights(0); mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast(iteration) + "_" + name + ".fib"); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); peak_image_writer->Update(); std::cout.rdbuf (old); // <-- restore // logfile << name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; } } clock.Stop(); int h = static_cast(clock.GetTotal()/3600); int m = (static_cast(clock.GetTotal())%3600)/60; int s = static_cast(clock.GetTotal())%60; MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s"; logfile.close(); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/CalculateOverlap.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CalculateOverlap.cpp index 1353955123..1c223ebadd 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CalculateOverlap.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CalculateOverlap.cpp @@ -1,108 +1,108 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUIntImgType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Calculate Overlap"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("tractogram", "", mitkCommandLineParser::String, "", "", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("mask", "", mitkCommandLineParser::String, "", "", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("peaks", "", mitkCommandLineParser::String, "", "", us::Any(), true, false, false, mitkCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string input_tractogram = us::any_cast(parsedArgs["tractogram"]); std::string mask_file = us::any_cast(parsedArgs["mask"]); try { mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(input_tractogram); mitk::Image::Pointer mask = mitk::IOUtil::Load(mask_file); ItkFloatImgType::Pointer itk_mask; mitk::CastToItkImage(mask, itk_mask); mitk::PeakImage::ItkPeakImageType::Pointer peaks = nullptr; if (parsedArgs.count("peaks")) { mitk::Image::Pointer mitk_peaks = mitk::IOUtil::Load(us::any_cast(parsedArgs["peaks"])); typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitk_peaks); caster->Update(); peaks = caster->GetOutput(); float overlap = 0; float directional_overlap = 0; std::tie(directional_overlap, overlap) = fib->GetDirectionalOverlap(itk_mask, peaks); MITK_INFO << "Overlap<<" << overlap; MITK_INFO << "DirectionalOverlap<<" << directional_overlap; } else MITK_INFO << "Overlap<<" << fib->GetOverlap(itk_mask); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/CheckEpsAndOverlap.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CheckEpsAndOverlap.cpp index 491d3e675f..8f0bb5d803 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CheckEpsAndOverlap.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CheckEpsAndOverlap.cpp @@ -1,98 +1,98 @@ /*=================================================================== 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 typedef itksys::SystemTools ist; typedef itk::Image ItkUcharImgType; typedef itk::Image ItkUIntImgType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Calculate Overlap"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("tractogram", "", mitkCommandLineParser::String, "", "file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("overlap_image", "", mitkCommandLineParser::String, "", "file", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("ep_image", "", mitkCommandLineParser::String, "", "file", us::Any(), false, false, false, mitkCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string input_tractogram = us::any_cast(parsedArgs["tractogram"]); std::string mask_file = us::any_cast(parsedArgs["overlap_image"]); std::string mask_file2 = us::any_cast(parsedArgs["ep_image"]); try { mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(input_tractogram); mitk::Image::Pointer mask = mitk::IOUtil::Load(mask_file); mitk::Image::Pointer mask2 = mitk::IOUtil::Load(mask_file2); ItkUcharImgType::Pointer itk_mask; mitk::CastToItkImage(mask, itk_mask); ItkUcharImgType::Pointer itk_mask2; mitk::CastToItkImage(mask2, itk_mask2); float ol = fib->GetOverlap(itk_mask); float ep = fib->GetNumEpFractionInMask(itk_mask2, true); MITK_INFO << "Overlap<<" << ol; MITK_INFO << "EP-Fraction<<" << ep; } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp index 6b7313bcff..d04dcb3773 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp @@ -1,199 +1,199 @@ /*=================================================================== 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 typedef itk::Image ItkFloatImgType; /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Extract Similar Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle (.fib, .trk, .tck)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("ref_tracts", "", mitkCommandLineParser::StringList, "Ref. Tracts:", "reference tracts (.fib, .trk, .tck)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("ref_masks", "", mitkCommandLineParser::StringList, "Ref. Masks:", "reference bundle masks", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("distance", "", mitkCommandLineParser::Int, "Distance:", "", 10); parser.addArgument("metric", "", mitkCommandLineParser::String, "Metric:", "EU_MEAN (default), EU_STD, EU_MAX"); parser.addArgument("subsample", "", mitkCommandLineParser::Float, "Subsampling factor:", "Only use specified fraction of input fibers", 1.0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string in_fib = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType ref_bundle_files = us::any_cast(parsedArgs["ref_tracts"]); mitkCommandLineParser::StringContainerType ref_mask_files; if (parsedArgs.count("ref_masks")) ref_mask_files = us::any_cast(parsedArgs["ref_masks"]); if (ref_mask_files.size()>0 && ref_mask_files.size()!=ref_bundle_files.size()) { MITK_INFO << "If reference masks are used, there has to be one mask per reference tract."; return EXIT_FAILURE; } int distance = 10; if (parsedArgs.count("distance")) distance = us::any_cast(parsedArgs["distance"]); std::string metric = "EU_MEAN"; if (parsedArgs.count("metric")) metric = us::any_cast(parsedArgs["metric"]); float subsample = 1.0; if (parsedArgs.count("subsample")) subsample = us::any_cast(parsedArgs["subsample"]); try { mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(in_fib); std::srand(0); if (subsample<1.0f) fib = fib->SubsampleFibers(subsample, true); mitk::FiberBundle::Pointer resampled_fib = fib->GetDeepCopy(); resampled_fib->ResampleToNumPoints(12); auto ref_fibs = mitk::DiffusionDataIOHelper::load_fibs(ref_bundle_files); auto ref_masks = mitk::DiffusionDataIOHelper::load_itk_images(ref_mask_files); std::vector< float > distances; distances.push_back(distance); mitk::FiberBundle::Pointer extracted = mitk::FiberBundle::New(nullptr); unsigned int c = 0; for (auto ref_fib : ref_fibs) { MITK_INFO << "Extracting " << ist::GetFilenameName(ref_bundle_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect try { itk::TractClusteringFilter::Pointer segmenter = itk::TractClusteringFilter::New(); // calculate centroids from reference bundle { itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances({10,20,30}); clusterer->SetTractogram(ref_fib); clusterer->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); clusterer->SetMergeDuplicateThreshold(0.0); clusterer->Update(); std::vector tracts = clusterer->GetOutCentroids(); ref_fib = mitk::FiberBundle::New(nullptr); ref_fib = ref_fib->AddBundles(tracts); mitk::IOUtil::Save(ref_fib, out_root + "centroids_" + ist::GetFilenameName(ref_bundle_files.at(c))); segmenter->SetInCentroids(ref_fib); } // segment tract if (cSetFilterMask(ref_masks.at(c)); segmenter->SetOverlapThreshold(0.8f); } segmenter->SetDistances(distances); segmenter->SetTractogram(resampled_fib); segmenter->SetMergeDuplicateThreshold(0.0); segmenter->SetDoResampling(false); if (metric=="EU_MEAN") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); else if (metric=="EU_STD") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); else if (metric=="EU_MAX") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMax()}); segmenter->Update(); std::vector< std::vector< unsigned int > > clusters = segmenter->GetOutFiberIndices(); if (clusters.size()>0) { vtkSmartPointer weights = vtkSmartPointer::New(); mitk::FiberBundle::Pointer result = mitk::FiberBundle::New(nullptr); std::vector< mitk::FiberBundle::Pointer > result_fibs; for (unsigned int cluster_index=0; cluster_indexGeneratePolyDataByIds(clusters.at(cluster_index), weights))); result = result->AddBundles(result_fibs); extracted = extracted->AddBundle(result); mitk::IOUtil::Save(result, out_root + "extracted_" + ist::GetFilenameName(ref_bundle_files.at(c))); fib = mitk::FiberBundle::New(fib->GeneratePolyDataByIds(clusters.back(), weights)); resampled_fib = mitk::FiberBundle::New(resampled_fib->GeneratePolyDataByIds(clusters.back(), weights)); } } catch(itk::ExceptionObject& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.GetDescription(); } catch(std::exception& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.what(); } std::cout.rdbuf (old); // <-- restore if (fib->GetNumFibers()==0) break; ++c; } MITK_INFO << "Extracted streamlines: " << extracted->GetNumFibers(); mitk::IOUtil::Save(extracted, out_root + "extracted_streamlines.trk"); MITK_INFO << "Residual streamlines: " << fib->GetNumFibers(); mitk::IOUtil::Save(fib, out_root + "residual_streamlines.trk"); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp index f66fc49670..67bad8e9b2 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp @@ -1,167 +1,167 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Get Overlapping Tracts"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Find tracts that overlap with the reference masks or tracts"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::StringList, "Input:", "input tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output Folder:", "move input tracts that do/don't overlap here", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("reference", "", mitkCommandLineParser::StringList, "Reference:", "reference tractograms or mask images", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "", 0.9); parser.addArgument("use_any_overlap", "", mitkCommandLineParser::Bool, "Use any overlap:", "Don't find maximum overlap but use first overlap larger threshold"); parser.addArgument("dont_save_tracts", "", mitkCommandLineParser::Bool, "Don't save tracts:", "if true, only text files documenting the overlaps are saved and no tract files are copied"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input = us::any_cast(parsedArgs["i"]); mitkCommandLineParser::StringContainerType reference = us::any_cast(parsedArgs["reference"]); std::string out_folder = us::any_cast(parsedArgs["o"]); bool use_any_overlap = false; if (parsedArgs.count("use_any_overlap")) use_any_overlap = us::any_cast(parsedArgs["use_any_overlap"]); bool dont_save_tracts = false; if (parsedArgs.count("dont_save_tracts")) dont_save_tracts = us::any_cast(parsedArgs["dont_save_tracts"]); float overlap_threshold = 0.9; if (parsedArgs.count("overlap_fraction")) overlap_threshold = us::any_cast(parsedArgs["overlap_fraction"]); try { MITK_INFO << "Loading references"; std::vector< std::string > reference_names; auto masks = mitk::DiffusionDataIOHelper::load_itk_images(reference, &reference_names); auto reference_fibs = mitk::DiffusionDataIOHelper::load_fibs(reference, &reference_names); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect itk::TractDensityImageFilter< ItkFloatImgType >::Pointer filter = itk::TractDensityImageFilter< ItkFloatImgType >::New(); filter->SetUpsamplingFactor(0.25); filter->SetBinaryOutput(true); for (auto fib : reference_fibs) { filter->SetFiberBundle(fib); filter->Update(); masks.push_back(filter->GetOutput()); } std::cout.rdbuf (old); // <-- restore MITK_INFO << "Loading input tractograms"; std::vector< std::string > input_names; auto input_fibs = mitk::DiffusionDataIOHelper::load_fibs(input, &input_names); MITK_INFO << "Finding overlaps"; ofstream logfile; logfile.open (out_folder + "Overlaps.txt"); ofstream logfile2; logfile2.open (out_folder + "AllOverlaps.txt"); boost::progress_display disp(input.size()); unsigned int c = 0; for (auto fib : input_fibs) { ++disp; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect bool is_overlapping = false; float overlap = 0; float max_overlap = 0; std::string max_ref = "-"; int i = 0; std::string overlap_string = ist::GetFilenameWithoutExtension(input_names.at(c)); for (auto m : masks) { overlap = fib->GetOverlap(m); if (overlap>max_overlap) { max_overlap = overlap; max_ref = ist::GetFilenameWithoutExtension(reference_names.at(i)); } if (use_any_overlap && overlap>=overlap_threshold) break; overlap_string += " " + ist::GetFilenameWithoutExtension(reference_names.at(i)) + " " + boost::lexical_cast(overlap); ++i; } if (overlap>=overlap_threshold) is_overlapping = true; logfile << ist::GetFilenameWithoutExtension(input_names.at(c)) << " - " << max_ref << ": " << boost::lexical_cast(max_overlap) << "\n"; logfile2 << overlap_string << "\n"; if (!dont_save_tracts && is_overlapping) ist::CopyAFile(input_names.at(c), out_folder + ist::GetFilenameName(input_names.at(c))); std::cout.rdbuf (old); // <-- restore ++c; } logfile.close(); logfile2.close(); } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp index c6981e34e8..93042bbdcc 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp @@ -1,214 +1,214 @@ /*=================================================================== 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 typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUIntImgType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Merge Overlapping Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::StringList, "Input:", "input tracts", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output Folder:", "output folder", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Tracts with overlap larger than this threshold are merged", 0.8, false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_folder = us::any_cast(parsedArgs["i"]); std::string out_folder = us::any_cast(parsedArgs["o"]); float overlap = 0.8; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); try { if (!ist::PathExists(out_folder)) ist::MakeDirectory(out_folder); std::vector< mitk::FiberBundle::Pointer > fibs = mitk::DiffusionDataIOHelper::load_fibs(input_folder); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::FiberBundle::Pointer combined = mitk::FiberBundle::New(); combined = combined->AddBundles(fibs); itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::Pointer endings = itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::New(); endings->SetFiberBundle(combined); endings->SetUpsamplingFactor(0.25); endings->Update(); ItkFloatImgType::Pointer ref_image = endings->GetOutput(); std::cout.rdbuf (old); // <-- restore for (int its = 0; its<3; its++) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect std::vector< ItkFloatImgType::Pointer > mask_images; for (auto fib : fibs) { itk::TractDensityImageFilter< ItkFloatImgType >::Pointer masks = itk::TractDensityImageFilter< ItkFloatImgType >::New(); masks->SetInputImage(ref_image); masks->SetBinaryOutput(true); masks->SetFiberBundle(fib); masks->SetUseImageGeometry(true); masks->Update(); mask_images.push_back(masks->GetOutput()); } int r=0; vnl_matrix< int > mat; mat.set_size(mask_images.size(), mask_images.size()); mat.fill(0); for (auto m1 : mask_images) { float max_overlap = overlap; int c = 0; for (auto m2 : mask_images) { if (c<=r) { ++c; continue; } itk::ImageRegionConstIterator it1(m1, m1->GetLargestPossibleRegion()); itk::ImageRegionConstIterator it2(m2, m2->GetLargestPossibleRegion()); unsigned int c1 = 0; unsigned int c2 = 0; unsigned int intersect = 0; while( !it1.IsAtEnd() ) { if( it1.Get()>0 && it2.Get()>0) ++intersect; if(it1.Get()>0) ++c1; if(it2.Get()>0) ++c2; ++it1; ++it2; } if ( (float)intersect/c1>max_overlap ) { max_overlap = (float)intersect/c1; mat.put(r,c, 1); } if ( (float)intersect/c2>max_overlap ) { max_overlap = (float)intersect/c2; mat.put(r,c, 1); } ++c; } ++r; } std::vector< mitk::FiberBundle::Pointer > out_fibs; std::vector< bool > used; for (unsigned int i=0; i0) { fib = fib->AddBundle(fibs.at(c)); used[c] = true; } } out_fibs.push_back(fib); } std::cout.rdbuf (old); // <-- restore MITK_INFO << fibs.size() << " --> " << out_fibs.size(); if (fibs.size()==out_fibs.size()) break; fibs = out_fibs; } int c = 0; for (auto fib : fibs) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + "/bundle_" + boost::lexical_cast(c) + ".trk"); std::cout.rdbuf (old); // <-- restore ++c; } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp index b5a8a7b3c3..57a3081893 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp @@ -1,190 +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 #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include typedef itk::Image< unsigned char, 3 > ItkUcharImageType; /*! \brief Calculate angular error between two sets of directions stored in multiple 3D vector images where each pixel corresponds to a vector (itk::Image< itk::Vector< float, 3>, 3 >) */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("test", "", mitkCommandLineParser::StringList, "Test images", "test direction images", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("reference", "", mitkCommandLineParser::StringList, "Reference images", "reference direction images", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output folder", "output folder", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("masks", "", mitkCommandLineParser::StringList, "Mask(s)", "mask image(s)", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("verbose", "", mitkCommandLineParser::Bool, "Verbose", "output error images"); parser.addArgument("ignore_test", "", mitkCommandLineParser::Bool, "Ignore missing test", "don't increase error if no test directions are found"); parser.addArgument("ignore_ref", "", mitkCommandLineParser::Bool, "Ignore ignore missing ref", "don't increase error if no ref directions are found"); parser.setCategory("Fiber Tracking Evaluation"); parser.setTitle("Peaks Angular Error"); parser.setDescription("Calculate angular error between two sets of peak images (1-1 correspondence)"); parser.setContributor("MIC"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType testImages = us::any_cast(parsedArgs["test"]); mitkCommandLineParser::StringContainerType referenceImages = us::any_cast(parsedArgs["reference"]); mitkCommandLineParser::StringContainerType maskImages; if (parsedArgs.count("masks")) maskImages = us::any_cast(parsedArgs["masks"]); std::string outRoot = us::any_cast(parsedArgs["o"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); bool ignore_test = false; if (parsedArgs.count("ignore_test")) ignore_test = us::any_cast(parsedArgs["ignore_test"]); bool ignore_ref = false; if (parsedArgs.count("ignore_ref")) ignore_ref = us::any_cast(parsedArgs["ignore_ref"]); try { typedef itk::ComparePeakImagesFilter< float > EvaluationFilterType; std::vector test_names; auto test_images = mitk::DiffusionDataIOHelper::load_itk_images(testImages, &test_names); // load reference directions std::vector ref_names; auto ref_images = mitk::DiffusionDataIOHelper::load_itk_images(referenceImages, &ref_names); // load/create mask image auto itkMaskImages = mitk::DiffusionDataIOHelper::load_itk_images(maskImages); if (test_images.size()!=ref_images.size()) mitkThrow() << "Matching number of test and reference image required!"; for (unsigned int i=0; iSetTestImage(test_images.at(i)); evaluationFilter->SetReferenceImage(ref_images.at(i)); if (iSetMaskImage(itkMaskImages.at(i)); evaluationFilter->SetIgnoreMissingTestDirections(ignore_test); evaluationFilter->SetIgnoreMissingRefDirections(ignore_ref); evaluationFilter->Update(); std::string ref_name = ist::GetFilenameWithoutExtension(ref_names.at(i)); std::string test_name = ist::GetFilenameWithoutExtension(test_names.at(i)); if (verbose) { mitk::LocaleSwitch localeSwitch("C"); EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); EvaluationFilterType::OutputImageType::Pointer lengthErrorImage = evaluationFilter->GetOutput(1); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; { WriterType::Pointer writer = WriterType::New(); std::string outfilename = outRoot; outfilename.append(ref_name + "_" + test_name + "_AngularError.nii.gz"); writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } { WriterType::Pointer writer = WriterType::New(); std::string outfilename = outRoot; outfilename.append(ref_name + "_" + test_name + "_LengthError.nii.gz"); writer->SetFileName(outfilename.c_str()); writer->SetInput(lengthErrorImage); writer->Update(); } } std::string logFile = outRoot; logFile.append("AngularErrors.csv"); bool add_header = true; if (ist::FileExists(logFile, true)) add_header = false; ofstream file; file.open (logFile.c_str(), std::fstream::app); std::string sens; if (add_header) sens.append("Test,Reference,Mean,Median,Maximum,Minimum,Stdev\n"); sens.append(test_name); sens.append(","); sens.append(ref_name); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append("\n"); std::cout << sens; file << sens; file.close(); } } catch (const itk::ExceptionObject& e) { - std::cout << e; + std::cout << e.what(); 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/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp index 4318ad6ac6..37e07f888e 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp @@ -1,149 +1,149 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include typedef itk::Image< unsigned char, 3 > ItkUcharImageType; typedef itk::Image< float, 4 > ItkPeakImgType; int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Reference Similarity"); parser.setCategory("Fiber Tracking Evaluation"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::StringList, "Input Tracts:", "input tracts folder", us::Any(), false); parser.addArgument("reference_tracts", "", mitkCommandLineParser::StringList, "", "", us::Any(), false); parser.addArgument("reference_masks", "", mitkCommandLineParser::StringList, "", "", us::Any(), false); parser.addArgument("reference_peaks", "", mitkCommandLineParser::StringList, "", "", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "", "", us::Any(), false); parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 20); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string out_folder = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType input_tract_files = us::any_cast(parsedArgs["i"]); mitkCommandLineParser::StringContainerType reference_tract_files = us::any_cast(parsedArgs["reference_tracts"]); mitkCommandLineParser::StringContainerType reference_mask_files = us::any_cast(parsedArgs["reference_masks"]); mitkCommandLineParser::StringContainerType reference_peak_files = us::any_cast(parsedArgs["reference_peaks"]); int fiber_points = 20; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); try { std::vector input_tract_names; std::vector ref_tract_names; std::vector< mitk::FiberBundle::Pointer > input_tracts = mitk::DiffusionDataIOHelper::load_fibs(input_tract_files, &input_tract_names); std::vector< mitk::FiberBundle::Pointer > reference_tracts = mitk::DiffusionDataIOHelper::load_fibs(reference_tract_files, &ref_tract_names); std::vector< ItkUcharImageType::Pointer > reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files); std::vector< ItkPeakImgType::Pointer > reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peak_files); MITK_INFO << "Calculating distances"; itk::TractDistanceFilter::Pointer distance_calculator = itk::TractDistanceFilter::New(); distance_calculator->SetNumPoints(fiber_points); distance_calculator->SetTracts1(input_tracts); distance_calculator->SetTracts2(reference_tracts); distance_calculator->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); distance_calculator->Update(); auto distances = distance_calculator->GetAllDistances(); vnl_matrix voxel_overlap; voxel_overlap.set_size(input_tracts.size(), reference_tracts.size()); vnl_matrix dir_overlap; dir_overlap.set_size(input_tracts.size(), reference_tracts.size()); MITK_INFO << "Calculating overlap"; boost::progress_display disp(input_tracts.size()*reference_tracts.size()); int r=0; for (auto fib : input_tracts) { int c=0; for (auto ref_mask : reference_masks) { // ++disp; // std::streambuf *old = cout.rdbuf(); // <-- save // std::stringstream ss; // std::cout.rdbuf (ss.rdbuf()); // <-- redirect float overlap = 0; float directional_overlap = 0; std::tie(directional_overlap, overlap) = fib->GetDirectionalOverlap(ref_mask, reference_peaks.at(c)); voxel_overlap[r][c] = overlap; dir_overlap[r][c] = directional_overlap; // std::cout.rdbuf (old); // <-- restore ++c; } ++r; } ofstream logfile; logfile.open(out_folder + "ref_tract_names.txt"); for (unsigned int i=0; i #include #include #include #include #include #include #include #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Tract Distance"); parser.setCategory("Fiber Processing"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkCommandLineParser::StringList, "Input tracts 1:", "input tracts 1", us::Any(), false); parser.addArgument("", "i2", mitkCommandLineParser::StringList, "Input tracts 2:", "input tracts 2", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output logfile", us::Any(), false); parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 12); parser.addArgument("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN (default), EU_STD, EU_MAX"); parser.addArgument("metric_weights", "", mitkCommandLineParser::StringList, "Metric weights:", "add one float weight for each used metric"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType t1_folder = us::any_cast(parsedArgs["i1"]); mitkCommandLineParser::StringContainerType t2_folder = us::any_cast(parsedArgs["i2"]); std::string out_file = us::any_cast(parsedArgs["o"]); int fiber_points = 12; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); std::vector< std::string > metric_strings = {"EU_MEAN"}; if (parsedArgs.count("metrics")) metric_strings = us::any_cast(parsedArgs["metrics"]); std::vector< std::string > metric_weights = {"1.0"}; if (parsedArgs.count("metric_weights")) metric_weights = us::any_cast(parsedArgs["metric_weights"]); if (metric_strings.size()!=metric_weights.size()) { MITK_INFO << "Each metric needs an associated metric weight!"; return EXIT_FAILURE; } try { std::vector t1_files; std::vector< mitk::FiberBundle::Pointer > tractograms1 = mitk::DiffusionDataIOHelper::load_fibs(t1_folder, &t1_files); std::vector t2_files; std::vector< mitk::FiberBundle::Pointer > tractograms2 = mitk::DiffusionDataIOHelper::load_fibs(t2_folder, &t2_files); MITK_INFO << "Loaded " << tractograms1.size() << " source tractograms."; MITK_INFO << "Loaded " << tractograms2.size() << " target tractograms."; itk::TractDistanceFilter::Pointer distance_calculator = itk::TractDistanceFilter::New(); distance_calculator->SetNumPoints(fiber_points); distance_calculator->SetTracts1(tractograms1); distance_calculator->SetTracts2(tractograms2); std::vector< mitk::ClusteringMetric* > metrics; int mc = 0; for (auto m : metric_strings) { float w = boost::lexical_cast(metric_weights.at(mc)); MITK_INFO << "Metric: " << m << " (w=" << w << ")"; if (m=="EU_MEAN") metrics.push_back({new mitk::ClusteringMetricEuclideanMean()}); else if (m=="EU_STD") metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); else if (m=="EU_MAX") metrics.push_back({new mitk::ClusteringMetricEuclideanMax()}); metrics.back()->SetScale(w); mc++; } if (metrics.empty()) { MITK_INFO << "No metric selected!"; return EXIT_FAILURE; } distance_calculator->SetMetrics(metrics); distance_calculator->Update(); MITK_INFO << "Distances:"; auto distances = distance_calculator->GetMinDistances(); auto indices = distance_calculator->GetMinIndices(); ofstream logfile; logfile.open (out_file); for (unsigned int i=0; i