diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt index 6f1c0e87e9..6f29af6080 100755 --- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt @@ -1,67 +1,69 @@ OPTION(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF) IF(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS) # include necessary modules here MITK_CHECK_MODULE(_RESULT DiffusionCore FiberTracking ) IF(_RESULT) MESSAGE("Warning: DiffusionMiniApps is missing ${_RESULT}") ELSE(_RESULT) MITK_USE_MODULE( DiffusionCore FiberTracking ) # needed include directories INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${ALL_INCLUDE_DIRECTORIES}) PROJECT( mitkDiffusionMiniApps ) # fill in the standalone executables here SET(DIFFUSIONMINIAPPS mitkDiffusionMiniApps ) # set additional files here SET(DIFFUSIONCORE_ADDITIONAL_FILES MiniAppManager.cpp FileFormatConverter.cpp TensorReconstruction.cpp QballReconstruction.cpp DiffusionIndices.cpp CopyGeometry.cpp GibbsTracking.cpp StreamlineTracking.cpp FiberProcessing.cpp LocalDirectionalFiberPlausibility.cpp #TractogramAngularError.cpp FiberDirectionExtraction.cpp PeakExtraction.cpp PeaksAngularError.cpp MultishellMethods.cpp FiberFoxProcessing.cpp ExportShImage.cpp + NetworkCreation.cpp + NetworkStatistics.cpp ) # deprecated # FOREACH(tool ${DIFFUSIONMINIAPPS}) # ADD_EXECUTABLE( # ${tool} # ${tool}.cpp # ${DIFFUSIONCORE_ADDITIONAL_FILES} # ) # TARGET_LINK_LIBRARIES( # ${tool} # ${ALL_LIBRARIES} ) # ENDFOREACH(tool) mitk_create_executable(mitkDiffusionMiniApps - DEPENDS DiffusionCore FiberTracking + DEPENDS DiffusionCore FiberTracking Connectomics ) ENDIF() MITK_INSTALL_TARGETS(EXECUTABLES mitkDiffusionMiniApps ) ENDIF(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS) diff --git a/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp b/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp new file mode 100644 index 0000000000..e7b488d94d --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/NetworkCreation.cpp @@ -0,0 +1,145 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" + +// std includes +#include + +// CTK includes +#include "ctkCommandLineParser.h" + +// MITK includes +#include +#include +#include +#include +#include "mitkConnectomicsNetworkCreator.h" + +int NetworkCreation(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("fiberImage", "f", ctkCommandLineParser::String, "input fiber image (.fib)", us::Any(), false); + parser.addArgument("parcellation", "p", ctkCommandLineParser::String, "parcellation image", us::Any(), false); + parser.addArgument("outputNetwork", "o", ctkCommandLineParser::String, "where to save the ouput (.cnf)", us::Any(), false); + + parser.addArgument("radius", "r", ctkCommandLineParser::Int, "Search radius in mm", 15, true); + parser.addArgument("noCenterOfMass", "com", ctkCommandLineParser::Bool, "Do not use center of mass for node positions"); + + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + //default values + int searchRadius( 15 ); + bool noCenterOfMass( false ); + + // parse command line arguments + std::string fiberFilename = us::any_cast(parsedArgs["fiberImage"]); + std::string parcellationFilename = us::any_cast(parsedArgs["parcellation"]); + std::string outputFilename = us::any_cast(parsedArgs["outputNetwork"]); + + if (parsedArgs.count("radius")) + searchRadius = us::any_cast(parsedArgs["radius"]); + + + if (parsedArgs.count("noCenterOfMass")) + noCenterOfMass = us::any_cast(parsedArgs["noCenterOfMass"]); + + try + { + // registering Factories + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + RegisterConnectomicsObjectFactory(); + + const std::string s1="", s2=""; + + // load fiber image + std::vector fiberInfile = + mitk::BaseDataIO::LoadBaseDataFromFile( fiberFilename, s1, s2, false ); + if( fiberInfile.empty() ) + { + std::string errorMessage = "Fiber Image at " + fiberFilename + " could not be read. Aborting."; + MITK_ERROR << errorMessage; + return EXIT_FAILURE; + } + mitk::BaseData* fiberBaseData = fiberInfile.at(0); + mitk::FiberBundleX* fiberBundle = dynamic_cast( fiberBaseData ); + + // load parcellation + std::vector parcellationInFile = + mitk::BaseDataIO::LoadBaseDataFromFile( parcellationFilename, s1, s2, false ); + if( parcellationInFile.empty() ) + { + std::string errorMessage = "Parcellation at " + parcellationFilename + " could not be read. Aborting."; + MITK_ERROR << errorMessage; + return EXIT_FAILURE; + } + mitk::BaseData* parcellationBaseData = parcellationInFile.at(0); + mitk::Image* parcellationImage = dynamic_cast( parcellationBaseData ); + + + + // do creation + mitk::ConnectomicsNetworkCreator::Pointer connectomicsNetworkCreator = mitk::ConnectomicsNetworkCreator::New(); + connectomicsNetworkCreator->SetSegmentation( parcellationImage ); + connectomicsNetworkCreator->SetFiberBundle( fiberBundle ); + if( !noCenterOfMass ) + { + connectomicsNetworkCreator->CalculateCenterOfMass(); + } + connectomicsNetworkCreator->SetEndPointSearchRadius( searchRadius ); + connectomicsNetworkCreator->CreateNetworkFromFibersAndSegmentation(); + + + mitk::ConnectomicsNetwork::Pointer network = connectomicsNetworkCreator->GetNetwork(); + + MITK_INFO << "searching writer"; + mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); + for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) + { + if ( (*it)->CanWriteBaseDataType(network.GetPointer()) ) + { + MITK_INFO << "writing"; + (*it)->SetFileName( outputFilename.c_str() ); + (*it)->DoWrite( network.GetPointer() ); + } + } + + return EXIT_SUCCESS; + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + MITK_INFO << "DONE"; + return EXIT_SUCCESS; +} +RegisterDiffusionMiniApp(NetworkCreation); diff --git a/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp new file mode 100644 index 0000000000..b9b9321541 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/NetworkStatistics.cpp @@ -0,0 +1,338 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" + +// std includes +#include +#include +#include +#include + +// ITK includes +#include + +// CTK includes +#include "ctkCommandLineParser.h" + +// MITK includes +#include +#include +#include +#include + +int NetworkStatistics(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("inputNetwork", "i", ctkCommandLineParser::String, "input connectomics network (.cnf)", us::Any(), false); + parser.addArgument("outputFile", "o", ctkCommandLineParser::String, "name of output file", us::Any(), false); + + parser.addArgument("noGlobalStatistics", "g", ctkCommandLineParser::Bool, "Do not calculate global statistics"); + parser.addArgument("createConnectivityMatriximage", "I", ctkCommandLineParser::Bool, "Write connectivity matrix image"); + parser.addArgument("binaryConnectivity", "b", ctkCommandLineParser::Bool, "Whether to create a binary connectivity matrix"); + parser.addArgument("rescaleConnectivity", "r", ctkCommandLineParser::Bool, "Whether to rescale the connectivity matrix"); + parser.addArgument("localStatistics", "L", ctkCommandLineParser::StringList, "Provide a list of node labels for local statistics", us::Any()); + + + 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 ); + + + // parse command line arguments + std::string networkName = us::any_cast(parsedArgs["inputNetwork"]); + std::string outName = us::any_cast(parsedArgs["outputFile"]); + + ctkCommandLineParser::StringContainerType localLabels; + + if(parsedArgs.count("localStatistics")) + { + localLabels = us::any_cast(parsedArgs["localStatistics"]); + } + + 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"]); + + try + { + const std::string s1="", s2=""; + RegisterConnectomicsObjectFactory(); + + // load network + std::vector networkFile = + mitk::BaseDataIO::LoadBaseDataFromFile( networkName, s1, s2, false ); + if( networkFile.empty() ) + { + std::string errorMessage = "File at " + networkName + " could not be read. Aborting."; + MITK_ERROR << errorMessage; + return EXIT_FAILURE; + } + mitk::BaseData* networkBaseData = networkFile.at(0); + mitk::ConnectomicsNetwork* network = dynamic_cast( networkBaseData ); + + if( !network ) + { + std::string errorMessage = "Read file at " + networkName + " could not be recognized as network. Aborting."; + MITK_ERROR << errorMessage; + return EXIT_FAILURE; + } + + mitk::ConnectomicsStatisticsCalculator::Pointer statisticsCalculator = mitk::ConnectomicsStatisticsCalculator::New(); + + statisticsCalculator->SetNetwork( network ); + statisticsCalculator->Update(); + + // global statistics + if( !noGlobalStatistics ) + { + std::string globalOutName = outName + "_global.txt"; + std::stringstream headerStream; + headerStream << "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 dataStream; + dataStream << 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; + + ofstream outFile( globalOutName.c_str(), ios::out ); + + if( ! outFile.is_open() ) + { + std::string errorMessage = "Could not open " + globalOutName + " for writing."; + MITK_ERROR << errorMessage; + return EXIT_FAILURE; + } + + outFile << headerStream.str() << dataStream.str(); + outFile.close(); + } // end global statistics + + //create connectivity matrix png + if( createConnectivityMatriximage ) + { + std::string connectivity_png_postfix = "_connectivity"; + if( binaryConnectivity ) + { + connectivity_png_postfix += "_binary"; + } + else if( rescaleConnectivity ) + { + connectivity_png_postfix += "_rescaled"; + } + connectivity_png_postfix += ".png"; + + /* File format + * A png file depicting the binary connectivity matrix + */ + itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::Pointer filter = itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::New(); + filter->SetInputNetwork( network ); + filter->SetBinaryConnectivity( binaryConnectivity ); + filter->SetRescaleConnectivity( rescaleConnectivity ); + filter->Update(); + + typedef itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::OutputImageType connectivityMatrixImageType; + + itk::ImageFileWriter< connectivityMatrixImageType >::Pointer connectivityWriter = itk::ImageFileWriter< connectivityMatrixImageType >::New(); + + connectivityWriter->SetInput( filter->GetOutput() ); + connectivityWriter->SetFileName( outName + connectivity_png_postfix); + connectivityWriter->Update(); + + MITK_INFO << "Connectivity matrix image written."; + } // end create connectivity matrix png + + // calculate local averages + if( localLabels.size() > 0 ) + { + std::string localOutName = outName + "_local.txt"; + + // Create LabelToIndex translation + std::map< std::string, int > labelToIdMap; + std::vector< mitk::ConnectomicsNetwork::NetworkNode > nodeVector = network->GetVectorOfAllNodes(); + for(int loop(0); loop < nodeVector.size(); loop++) + { + labelToIdMap.insert( std::pair< std::string, int>(nodeVector.at(loop).label, nodeVector.at(loop).id) ); + } + + std::vector< int > degreeVector = network->GetDegreeOfNodes(); + std::vector< double > ccVector = network->GetLocalClusteringCoefficients( ); + std::vector< double > bcVector = network->GetNodeBetweennessVector( ); + + double sumDegree( 0 ); + double sumCC( 0 ); + double sumBC( 0 ); + double count( 0 ); + + for( int loop(0); loop < localLabels.size(); loop++ ) + { + if( network->CheckForLabel(localLabels.at( loop )) ) + { + sumDegree = sumDegree + degreeVector.at( labelToIdMap.find( localLabels.at( loop ) )->second ); + sumCC = sumCC + ccVector.at( labelToIdMap.find( localLabels.at( loop ) )->second ); + sumBC = sumBC + bcVector.at( labelToIdMap.find( localLabels.at( loop ) )->second ); + count = count + 1; + } + else + { + MITK_ERROR << "Illegal label. Label: \"" << localLabels.at( loop ) << "\" not found."; + } + } + + std::stringstream headerStream; + headerStream << "LocalAverageDegree " + << "LocalAverageCC " + << "LocalAverageBC " + << "NumberOfNodes" + << std::endl; + + std::stringstream dataStream; + dataStream << sumDegree / count << " " + << sumCC / count << " " + << sumBC / count << " " + << count + << std::endl; + + ofstream outFile( localOutName.c_str(), ios::out ); + + if( ! outFile.is_open() ) + { + std::string errorMessage = "Could not open " + localOutName + " for writing."; + MITK_ERROR << errorMessage; + return EXIT_FAILURE; + } + + outFile << headerStream.str() << dataStream.str(); + outFile.close(); + }// end calculate local averages + + return EXIT_SUCCESS; + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + MITK_INFO << "DONE"; + return EXIT_SUCCESS; +} +RegisterDiffusionMiniApp(NetworkStatistics); diff --git a/Modules/DiffusionImaging/MiniApps/files.cmake b/Modules/DiffusionImaging/MiniApps/files.cmake index 45b39388e7..bebe55de6e 100644 --- a/Modules/DiffusionImaging/MiniApps/files.cmake +++ b/Modules/DiffusionImaging/MiniApps/files.cmake @@ -1,20 +1,22 @@ set(CPP_FILES mitkDiffusionMiniApps.cpp MiniAppManager.cpp FileFormatConverter.cpp TensorReconstruction.cpp QballReconstruction.cpp DiffusionIndices.cpp CopyGeometry.cpp GibbsTracking.cpp StreamlineTracking.cpp FiberProcessing.cpp LocalDirectionalFiberPlausibility.cpp #TractogramAngularError.cpp FiberDirectionExtraction.cpp PeakExtraction.cpp PeaksAngularError.cpp MultishellMethods.cpp FiberFoxProcessing.cpp ExportShImage.cpp + NetworkCreation.cpp + NetworkStatistics.cpp )