#include #include #include #include #include #include #include #include #include //////////////////////////////////////////////////////////////////////////////// // class NormalizedRootMeanSquareDifferenceCommand //////////////////////////////////////////////////////////////////////////////// gkg::NormalizedRootMeanSquareDifferenceCommand::NormalizedRootMeanSquareDifferenceCommand( int32_t argc, char* argv[], bool loadPlugin, bool removeFirst ) : gkg::Command( argc, argv, loadPlugin, removeFirst ) { try { parse(); } GKG_CATCH( "gkg::NormalizedRootMeanSquareDifferenceCommand::" "NormalizedRootMeanSquareDifferenceCommand( int32_t argc, " "char* argv[], bool loadPlugin, bool removeFirst ) " " : gkg::Command( argc, argv, loadPlugin, removeFirst )" ); } gkg::NormalizedRootMeanSquareDifferenceCommand::NormalizedRootMeanSquareDifferenceCommand( const std::vector< std::string >& fileNameDws, const std::string& fileNameDwMask, const std::vector< std::string >& fileNameT1s, const std::vector< std::string >& fileNameT2s, const std::string& fileNameRelaxometryMask, const std::string& fileNameOut, bool verbose ) { try { execute( fileNameDws, fileNameDwMask, fileNameT1s, fileNameT2s, fileNameRelaxometryMask, fileNameOut, verbose ); } GKG_CATCH( "gkg::NormalizedRootMeanSquareDifferenceCommand::" "NormalizedRootMeanSquareDifferenceCommand( " "const std::vector< std::string >& fileNameDws," "const std::string& fileNameDwMask," "const std::vector< std::string >& fileNameT1s," "const std::vector< std::string >& fileNameT2s," "const std::string& fileNameRelaxometryMask," "const std::string& fileNameOut," "bool verbose )" ); } gkg::NormalizedRootMeanSquareDifferenceCommand::NormalizedRootMeanSquareDifferenceCommand( const gkg::Dictionary& parameters ) : gkg::Command( parameters ) { try { DECLARE_VECTOR_OF_STRINGS_PARAMETER( parameters, std::vector< std::string >, fileNameDws ); DECLARE_STRING_PARAMETER( parameters, std::string, fileNameDwMask ); DECLARE_VECTOR_OF_STRINGS_PARAMETER( parameters, std::vector< std::string >, fileNameT1s ); DECLARE_VECTOR_OF_STRINGS_PARAMETER( parameters, std::vector< std::string >, fileNameT2s ); DECLARE_STRING_PARAMETER( parameters, std::string, fileNameRelaxometryMask ); DECLARE_STRING_PARAMETER( parameters, std::string, fileNameOut ); DECLARE_BOOLEAN_PARAMETER( parameters, bool, verbose ); execute( fileNameDws, fileNameDwMask, fileNameT1s, fileNameT2s, fileNameRelaxometryMask, fileNameOut, verbose ); } GKG_CATCH( "gkg::NormalizedRootMeanSquareDifferenceCommand::" "NormalizedRootMeanSquareDifferenceCommand( const " "gkg::Dictionary& parameters ) : gkg::Command( parameters )" ); } gkg::NormalizedRootMeanSquareDifferenceCommand:: ~NormalizedRootMeanSquareDifferenceCommand() { } std::string gkg::NormalizedRootMeanSquareDifferenceCommand::getStaticName() { try { return "NormalizedRootMeanSquareDifference"; } GKG_CATCH( "std::string gkg::" "NormalizedRootMeanSquareDifferenceCommand::getStaticName()" ); } void gkg::NormalizedRootMeanSquareDifferenceCommand::parse() { try { std::vector< std::string > fileNameDws; std::string fileNameDwMask; std::vector< std::string > fileNameT1s; std::vector< std::string > fileNameT2s; std::string fileNameRelaxometryMask; std::string fileNameOut; bool verbose = false; gkg::Application application( _argc, _argv, "Normalized Root Mean Square Difference", _loadPlugin ); application.addSeriesOption( "-dw", "Diffusion-weighted volume file name list", fileNameDws, 2 ); application.addSingleOption( "-m_dw", "Mask volume file name for diffusion", fileNameDwMask ); application.addSeriesOption( "-t1", "T1 VFA-SPGR volume file name list", fileNameT1s, 2 ); application.addSeriesOption( "-t2", "T2 MSME volume file name list", fileNameT2s, 2 ); application.addSingleOption( "-m_relaxo", "Mask volume file name for relaxometry", fileNameRelaxometryMask ); application.addSingleOption( "-o", "Output text file name (.txt)", fileNameOut ); application.addSingleOption( "-verbose", "Show as much information as possible", verbose, true ); application.initialize(); execute( fileNameDws, fileNameDwMask, fileNameT1s, fileNameT2s, fileNameRelaxometryMask, fileNameOut, verbose ); } GKG_CATCH_FUNCTOR( "void gkg::" "NormalizedRootMeanSquareDifferenceCommand::parse()" ); } void gkg::NormalizedRootMeanSquareDifferenceCommand::execute( const std::vector< std::string >& fileNameDws, const std::string& fileNameDwMask, const std::vector< std::string >& fileNameT1s, const std::vector< std::string >& fileNameT2s, const std::string& fileNameRelaxometryMask, const std::string& fileNameOut, bool verbose ) { try { //////////////////////////////////////////////////////////////////////////// // reading input DW files //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "reading '" << fileNameDws[ 0 ] << "' : " << std::flush; } gkg::Volume< float > firstDwVolumeIn; gkg::Reader::getInstance().read( fileNameDws[ 0 ], firstDwVolumeIn ); if ( verbose ) { std::cout << "done" << std::endl; } if ( verbose ) { std::cout << "reading '" << fileNameDws[ 1 ] << "' : " << std::flush; } gkg::Volume< float > secondDwVolumeIn; gkg::Reader::getInstance().read( fileNameDws[ 1 ], secondDwVolumeIn ); if ( verbose ) { std::cout << "done" << std::endl; } if ( !firstDwVolumeIn.hasSameSize( secondDwVolumeIn ) || !firstDwVolumeIn.hasSameResolution( secondDwVolumeIn, true ) ) { throw std::runtime_error( "incoherent size and resolution between input DW volume(s)" ); } //////////////////////////////////////////////////////////////////////////// // reading DW mask //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "reading '" << fileNameDwMask << "' : " << std::flush; } gkg::Volume< int16_t > maskDwVolume; gkg::Reader::getInstance().read( fileNameDwMask, maskDwVolume ); if ( verbose ) { std::cout << "done" << std::endl; } //////////////////////////////////////////////////////////////////////////// // performing normalized root mean square difference //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "Performing normalized RMS difference for DW data" << std::endl; } float dwSumOfSquare = 0.0f; float dwSumOfSquareDifference = 0.0f; float dwNormalizedRMSDifference = 0.0f; int32_t dwSizeX = firstDwVolumeIn.getSizeX(); int32_t dwSizeY = firstDwVolumeIn.getSizeY(); int32_t dwSizeZ = firstDwVolumeIn.getSizeZ(); int32_t dwSizeT = firstDwVolumeIn.getSizeT(); int32_t x = 0; int32_t y = 0; int32_t z = 0; int32_t t = 0; for ( z = 0; z < dwSizeZ; z++ ) { for ( y = 0; y < dwSizeY; y++ ) { for ( x = 0; x < dwSizeX; x++ ) { if ( maskDwVolume( x, y, z ) == 1 ) { for ( t = 0; t < dwSizeT; t++ ) { dwSumOfSquare += ( firstDwVolumeIn( x, y, z, t ) + \ secondDwVolumeIn( x, y, z, t ) / 2.0 ) * \ ( firstDwVolumeIn( x, y, z, t ) + \ secondDwVolumeIn( x, y, z, t ) / 2.0 ); dwSumOfSquareDifference += ( firstDwVolumeIn( x, y, z, t ) - \ secondDwVolumeIn( x, y, z, t ) ) * \ ( firstDwVolumeIn( x, y, z, t ) - \ secondDwVolumeIn( x, y, z, t ) ); } } } } } // std::cout << "sumOfSquare : " << dwSumOfSquare << std::endl; // std::cout << "sumOfSquareDifference : " // << dwSumOfSquareDifference // << std::endl; if ( dwSumOfSquare != 0 ) { dwNormalizedRMSDifference = sqrt( dwSumOfSquareDifference / dwSumOfSquare ); } std::cout << "===> normalizedRMSDifference : " << dwNormalizedRMSDifference << std::endl; //////////////////////////////////////////////////////////////////////////// // reading input T1 files //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "reading '" << fileNameT1s[ 0 ] << "' : " << std::flush; } gkg::Volume< float > firstT1VolumeIn; gkg::Reader::getInstance().read( fileNameT1s[ 0 ], firstT1VolumeIn ); if ( verbose ) { std::cout << "done" << std::endl; } if ( verbose ) { std::cout << "reading '" << fileNameT1s[ 1 ] << "' : " << std::flush; } gkg::Volume< float > secondT1VolumeIn; gkg::Reader::getInstance().read( fileNameT1s[ 1 ], secondT1VolumeIn ); if ( verbose ) { std::cout << "done" << std::endl; } if ( !firstT1VolumeIn.hasSameSize( secondT1VolumeIn ) || !firstT1VolumeIn.hasSameResolution( secondT1VolumeIn, true ) ) { throw std::runtime_error( "incoherent size and resolution between input T1 volume(s)" ); } //////////////////////////////////////////////////////////////////////////// // reading input T2 files //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "reading '" << fileNameT2s[ 0 ] << "' : " << std::flush; } gkg::Volume< float > firstT2VolumeIn; gkg::Reader::getInstance().read( fileNameT2s[ 0 ], firstT2VolumeIn ); if ( verbose ) { std::cout << "done" << std::endl; } if ( verbose ) { std::cout << "reading '" << fileNameT2s[ 1 ] << "' : " << std::flush; } gkg::Volume< float > secondT2VolumeIn; gkg::Reader::getInstance().read( fileNameT2s[ 1 ], secondT2VolumeIn ); if ( verbose ) { std::cout << "done" << std::endl; } if ( !firstT2VolumeIn.hasSameSize( secondT2VolumeIn ) || !firstT2VolumeIn.hasSameResolution( secondT2VolumeIn, true ) ) { throw std::runtime_error( "incoherent size and resolution between input T2 volume(s)" ); } //////////////////////////////////////////////////////////////////////////// // reading Relaxometry mask //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "reading '" << fileNameRelaxometryMask << "' : " << std::flush; } gkg::Volume< int16_t > maskRelaxometryVolume; gkg::Reader::getInstance().read( fileNameRelaxometryMask, maskRelaxometryVolume ); if ( verbose ) { std::cout << "done" << std::endl; } //////////////////////////////////////////////////////////////////////////// // performing normalized root mean square difference //////////////////////////////////////////////////////////////////////////// if ( verbose ) { std::cout << "Performing normalized RMS difference for T1 and T2 data" << std::endl; } float t1SumOfSquare = 0.0f; float t1SumOfSquareDifference = 0.0f; float t1NormalizedRMSDifference = 0.0f; float t2SumOfSquare = 0.0f; float t2SumOfSquareDifference = 0.0f; float t2NormalizedRMSDifference = 0.0f; int32_t t1SizeX = firstT1VolumeIn.getSizeX(); int32_t t1SizeY = firstT1VolumeIn.getSizeY(); int32_t t1SizeZ = firstT1VolumeIn.getSizeZ(); int32_t t1SizeT = firstT1VolumeIn.getSizeT(); for ( z = 0; z < t1SizeZ; z++ ) { for ( y = 0; y < t1SizeY; y++ ) { for ( x = 0; x < t1SizeX; x++ ) { if ( maskRelaxometryVolume( x, y, z ) == 1 ) { for ( t = 0; t < t1SizeT; t++ ) { t1SumOfSquare += ( firstT1VolumeIn( x, y, z, t ) + \ secondT1VolumeIn( x, y, z, t ) / 2.0 ) * \ ( firstT1VolumeIn( x, y, z, t ) + \ secondT1VolumeIn( x, y, z, t ) / 2.0 ); t1SumOfSquareDifference += ( firstT1VolumeIn( x, y, z, t ) - \ secondT1VolumeIn( x, y, z, t ) ) * \ ( firstT1VolumeIn( x, y, z, t ) - \ secondT1VolumeIn( x, y, z, t ) ); t2SumOfSquare += ( firstT2VolumeIn( x, y, z, t ) + \ secondT2VolumeIn( x, y, z, t ) / 2.0 ) * \ ( firstT2VolumeIn( x, y, z, t ) + \ secondT2VolumeIn( x, y, z, t ) / 2.0 ); t2SumOfSquareDifference += ( firstT2VolumeIn( x, y, z, t ) - \ secondT2VolumeIn( x, y, z, t ) ) * \ ( firstT2VolumeIn( x, y, z, t ) - \ secondT2VolumeIn( x, y, z, t ) ); } } } } } // std::cout << "T1 sumOfSquare : " << t1SumOfSquare << std::endl; // std::cout << "T1 sumOfSquareDifference : " // << t1SumOfSquareDifference // << std::endl; if ( t1SumOfSquare != 0 ) { t1NormalizedRMSDifference = sqrt( t1SumOfSquareDifference / t1SumOfSquare ); } std::cout << "===> T1 normalizedRMSDifference : " << t1NormalizedRMSDifference << std::endl; // std::cout << "T2 sumOfSquare : " << t2SumOfSquare << std::endl; // std::cout << "T2 sumOfSquareDifference : " // << t2SumOfSquareDifference // << std::endl; if ( t2SumOfSquare != 0 ) { t2NormalizedRMSDifference = sqrt( t2SumOfSquareDifference / t2SumOfSquare ); } std::cout << "===> T2 normalizedRMSDifference : " << t2NormalizedRMSDifference << std::endl; } GKG_CATCH( "void gkg::NormalizedRootMeanSquareDifferenceCommand::execute(" " const std::vector< std::string >& fileNameDws," " const std::string& fileNameDwMask," " const std::vector< std::string >& fileNameT1s," " const std::vector< std::string >& fileNameT2s," " const std::string& fileNameRelaxometryMask," " const std::string& fileNameOut," " bool verbose )" ); } RegisterCommandCreator( NormalizedRootMeanSquareDifferenceCommand, DECLARE_VECTOR_OF_STRINGS_PARAMETER_HELP( fileNameDws ) + DECLARE_STRING_PARAMETER_HELP( fileNameDwMask ) + DECLARE_VECTOR_OF_STRINGS_PARAMETER_HELP( fileNameT1s ) + DECLARE_VECTOR_OF_STRINGS_PARAMETER_HELP( fileNameT2s ) + DECLARE_STRING_PARAMETER_HELP( fileNameRelaxometryMask ) + DECLARE_STRING_PARAMETER_HELP( fileNameOut ) + DECLARE_BOOLEAN_PARAMETER_HELP( verbose ) );