diff --git a/Atlases.cmake b/Atlases.cmake new file mode 100644 index 0000000..b828ae6 --- /dev/null +++ b/Atlases.cmake @@ -0,0 +1,23 @@ +cmake_minimum_required(VERSION 2.7) + +include(${CMAKE_ROOT}/Modules/ExternalProject.cmake) + +set(ATLAS_URL "https://biomedic.doc.ic.ac.uk/brain-development/downloads/dHCP/atlases-dhcp-structural-pipeline-v1.zip") +set(ATLAS_MD5 77e924bc17a4906f5814874009f5eca6) + +if(NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/atlases) + ExternalProject_Add(atlases + URL ${ATLAS_URL} + URL_MD5 ${ATLAS_MD5} + PREFIX atlases + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + ) + + add_custom_target(atlases_move ALL + ${CMAKE_COMMAND} -E rename ${CMAKE_CURRENT_BINARY_DIR}/atlases/src/atlases ${CMAKE_CURRENT_SOURCE_DIR}/atlases + ) + + ADD_DEPENDENCIES(atlases_move atlases) +endif() diff --git a/CMakeLists.txt b/CMakeLists.txt index 8a8bcb9..e1a18e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,3 +25,6 @@ if (NOT COMMAND mirtk_configure_module) endif () mirtk_configure_module() + +SUBDIRS(ThirdParty/ANTs) +INCLUDE(Atlases.cmake) diff --git a/README.md b/README.md index 9086889..4ef27a1 100644 --- a/README.md +++ b/README.md @@ -1,58 +1,69 @@ -Draw-EM Segmentation Software -========================================== +# Draw-EM Segmentation Software + +![segmentation image](segmentation.png) Draw-EM (Developing brain Region Annotation With Expectation-Maximization) is a package of [MIRTK](https://github.com/BioMedIA/MIRTK) developed by Antonios Makropoulos and the [BioMedIA](https://biomedia.doc.ic.ac.uk/) research group. -It provides a collection of command-line tools as well as pipelines for the segmentation of developing brain MR images. +It provides a collection of command-line tools and pipelines for the segmentation of developing brain MR images. + +Draw-EM is used as part of the [dHCP structural pipeline](https://github.com/BioMedIA/dhcp-structural-pipeline) for the structural analysis (segmentation and surface extraction) of the neonatal brain. + + +## Dependencies +### FSL +The segmentation pipeline uses +[FSL](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL). +See the [installation instructions](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FslInstallation) for FSL. -Installation ------------- + +## Installation Draw-EM is part of MIRTK. -In order to compile it as part of MIRTK you need to do the following steps: -- download (clone) Draw-EM inside the Packages folder of your MIRTK directory -- enable compile of the package by setting the CMake flag "MODULE_DrawEM" of MIRTK to "ON" (using cmake or ccmake) +In order to compile it as part of MIRTK you need to: +- enable compilation of the package by setting the CMake flag "MODULE_DrawEM" of MIRTK to "ON" (using cmake or ccmake) - build MIRTK See the [installation instructions](https://mirtk.github.io/install.html) for a step-by-step guide on how to install MIRTK. -The segmentation pipeline further requires the following: -- [FSL](http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) installed -- The atlases required by Draw-EM need to be downloaded from [here](https://www.doc.ic.ac.uk/~am411/atlases-DrawEM.html) and extracted inside the Draw-EM directory. -- The N4 bias field correction from ITK is included in the ThirdParty/ITK folder. If the executable does not work, you will need to compile ITK and replace the ThirdParty/ITK/N4 binary with the N4BiasFieldCorrectionImageFilter binary - - -Run ---- +## Running the pipeline -The segmentation pipeline can be run with the following script: +The segmentation pipeline can be run as follows: -pipelines/neonatal-pipeline-v1.1.sh +mirtk neonatal-segmentation -The script requires the T2 image and the age at scan of the subject to be segmented (as first and second argument respectively). -Run the script without arguments for a detailed list of options. +``` +Arguments: + Nifti Image: The T2 image of the subject to be segmented. + Integer: Subject age in weeks. This is used to select the appropriate template for the initial registration. + If the age is <28w or >44w, it will be set to 28w or 44w respectively. +Options: + -d / -data-dir The directory used to run the script and output the files. + -c / -cleanup <0/1> Whether cleanup of temporary files is required (default: 1) + -p / -save-posteriors <0/1> Whether the structures' posteriors are required (default: 0) + -t / -threads Number of threads (CPU cores) allowed for the registration to run in parallel (default: 1) + -v / -verbose <0/1> Whether the script progress is reported (default: 1) + -h / -help / --help Print usage. +``` -License -------- +## License Draw-EM is distributed under the terms of the Apache License Version 2. See the accompanying [license file](LICENSE.txt) for details. The license enables usage of Draw-EM in both commercial and non-commercial applications, without restrictions on the licensing applied to the combined work. -Draw-EM uses third-party software, namely the "ITK: The Insight Toolkit for Segmentation and Registration". -ITK is distributed under the Apache License Version 2. -Specifically, the N4 bias field correction by Tustison et al. is included (http://www.insight-journal.org/browse/publication/640). -The covered file (N4) and license (LICENSE) can be found in ThirdParty/ITK. +## Releases +- v1.2: dHCP segmentation pipeline, method improvements described in [2]: multi-channel registration, modelling of hyper and hypo-intensities. +- v1.1: initial code release, method described in [1]. -Citation and acknowledgements ------------------------------ +## Citation and acknowledgements In case you found Draw-EM useful please give appropriate credit to the software. -Publication: +Publications: -A. Makropoulos et al. Automatic whole brain MRI segmentation of the developing neonatal brain, IEEE TMI, 2014 +1. A. Makropoulos et al. *"Automatic whole brain MRI segmentation of the developing neonatal brain"*, IEEE TMI, 2014 +2. A. Makropoulos, E. C. Robinson et al. *"The Developing Human Connectome Project: a Minimal Processing Pipeline for Neonatal Cortical Surface Reconstruction"*, NeuroImage, 2018 diff --git a/ThirdParty/ANTs/CMakeLists.txt b/ThirdParty/ANTs/CMakeLists.txt new file mode 100755 index 0000000..db43a6e --- /dev/null +++ b/ThirdParty/ANTs/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required(VERSION 2.7) +if(COMMAND cmake_policy) + cmake_policy(SET CMP0046 NEW) + cmake_policy(SET CMP0048 NEW) +endif(COMMAND cmake_policy) + +PROJECT( N4 ) + +IF(NOT ITK_DIR) + INCLUDE(External_ITK.cmake) + SET(External_ITK 1) +ENDIF() + +# Set up ITK +FIND_PACKAGE(ITK) +IF(ITK_FOUND) + INCLUDE(${ITK_USE_FILE}) +ELSE(ITK_FOUND) + MESSAGE(FATAL_ERROR + "Cannot build without ITK. Please set ITK_DIR.") +ENDIF(ITK_FOUND) + +INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS}) + +SET(SOURCES "${CMAKE_CURRENT_SOURCE_DIR}/antsCommandLineParser" "${CMAKE_CURRENT_SOURCE_DIR}/antsCommandLineOption" ) + +ADD_EXECUTABLE( N4 N4.cxx ${SOURCES} ) +TARGET_LINK_LIBRARIES(N4 ${ITK_LIBRARIES}) +set_target_properties(N4 PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/${MIRTK_TOOLS_DIR}") + +IF(External_ITK) + ADD_DEPENDENCIES(N4 ITK) +ENDIF() diff --git a/ThirdParty/ANTs/External_ITK.cmake b/ThirdParty/ANTs/External_ITK.cmake new file mode 100644 index 0000000..49d5fb5 --- /dev/null +++ b/ThirdParty/ANTs/External_ITK.cmake @@ -0,0 +1,18 @@ +cmake_minimum_required(VERSION 2.7) + +include(${CMAKE_ROOT}/Modules/ExternalProject.cmake) + +ExternalProject_Add( + ITK + GIT_REPOSITORY "https://github.com/InsightSoftwareConsortium/ITK.git" + SOURCE_DIR "${CMAKE_SOURCE_DIR}/ThirdParty/ITK" + CMAKE_ARGS -DBUILD_EXAMPLES=OFF -DBUILD_SHARED_LIBS=ON -DBUILD_TESTING=OFF -DCMAKE_INSTALL_PREFIX=${ITK_DIR} + UPDATE_COMMAND "" + PATCH_COMMAND "" + INSTALL_COMMAND "" +) + +ExternalProject_Get_Property(ITK SOURCE_DIR) +ExternalProject_Get_Property(ITK BINARY_DIR) +set(ITK_SOURCE_DIR ${SOURCE_DIR}) +set(ITK_DIR ${BINARY_DIR}) diff --git a/ThirdParty/ANTs/LICENSE.txt b/ThirdParty/ANTs/LICENSE.txt new file mode 100644 index 0000000..71e1383 --- /dev/null +++ b/ThirdParty/ANTs/LICENSE.txt @@ -0,0 +1,40 @@ +The file N4.cxx is based on the N4BiasFieldCorrection.cxx provided by ANTs. +Minor modifications have been applied to allow to build N4BiasFieldCorrection as an external package without the need to build the whole ANTs package. +Modifications are noted in the code with "modification{ ... }modification". + + +------------------------------------------------------------------------------ +The following copyright applies to the files: +antsCommandLineOption.cxx +antsCommandLineOption.h +antsCommandLineParser.cxx +antsCommandLineParser.h +N4BiasFieldCorrection.cxx + + +ConsortiumOfANTS® - http://www.picsl.upenn.edu/ANTS/ +Copyright (c) 2009-2013 (updated to current year ad infinitum) + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. +3. Neither the name of the consortium nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE CONSORTIUM AND CONTRIBUTORS ``AS IS'' AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +SUCH DAMAGE. diff --git a/ThirdParty/ANTs/N4.cxx b/ThirdParty/ANTs/N4.cxx new file mode 100755 index 0000000..a996c4e --- /dev/null +++ b/ThirdParty/ANTs/N4.cxx @@ -0,0 +1,1043 @@ +// This file is based on the N4BiasFieldCorrection.cxx provided by ANTs +// Minor modifications have been applied to allow to build N4BiasFieldCorrection as an external package withouth the need to build all of ANTs +// Modifications are noted in the following code with modification{ ... }modification + +//modification{ +/* +#include "antsUtilities.h" +#include "antsAllocImage.h" +*/ +//}modification +#include "antsCommandLineParser.h" + +//modification{ +//#include "ReadWriteData.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +//}modification + +#include "itkBinaryThresholdImageFilter.h" +#include "itkBSplineControlPointImageFilter.h" +#include "itkConstantPadImageFilter.h" +#include "itkExpImageFilter.h" +#include "itkExtractImageFilter.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkLabelStatisticsImageFilter.h" +#include "itkN4BiasFieldCorrectionImageFilter.h" +#include "itkShrinkImageFilter.h" +#include "itkTimeProbe.h" + +#include +#include +#include + +//modification{ +/* +#include "ANTsVersion.h" + +namespace ants +{ +*/ +//}modification +template +class CommandIterationUpdate : public itk::Command +{ +public: + typedef CommandIterationUpdate Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + itkNewMacro( Self ); +protected: + CommandIterationUpdate() + { + }; +public: + + void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE + { + Execute( (const itk::Object *) caller, event); + } + + void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE + { + const TFilter * filter = + dynamic_cast( object ); + + if( typeid( event ) != typeid( itk::IterationEvent ) ) + { + return; + } + if( filter->GetElapsedIterations() == 1 ) + { + std::cout << "Current level = " << filter->GetCurrentLevel() + 1 + << std::endl; + } + std::cout << " Iteration " << filter->GetElapsedIterations() + << " (of " + << filter->GetMaximumNumberOfIterations()[filter->GetCurrentLevel()] + << "). "; + std::cout << " Current convergence value = " + << filter->GetCurrentConvergenceMeasurement() + << " (threshold = " << filter->GetConvergenceThreshold() + << ")" << std::endl; + } +}; + +template +int N4( itk::ants::CommandLineParser *parser ) +{ + typedef float RealType; + + typedef itk::Image ImageType; + typename ImageType::Pointer inputImage = ITK_NULLPTR; + + typedef itk::Image MaskImageType; + typename MaskImageType::Pointer maskImage = ITK_NULLPTR; + + bool verbose = false; + typename itk::ants::CommandLineParser::OptionType::Pointer verboseOption = + parser->GetOption( "verbose" ); + if( verboseOption && verboseOption->GetNumberOfFunctions() ) + { + verbose = parser->Convert( verboseOption->GetFunction( 0 )->GetName() ); + } + + if( verbose ) + { + std::cout << std::endl << "Running N4 for " + << ImageDimension << "-dimensional images." << std::endl << std::endl; + } + + typedef itk::N4BiasFieldCorrectionImageFilter CorrecterType; + typename CorrecterType::Pointer correcter = CorrecterType::New(); + typename itk::ants::CommandLineParser::OptionType::Pointer inputImageOption = + parser->GetOption( "input-image" ); + +//modification{ + typedef itk::ImageFileReader ReaderType; + typename ReaderType::Pointer reader = ReaderType::New(); +//}modification + + if( inputImageOption && inputImageOption->GetNumberOfFunctions() ) + { + std::string inputFile = inputImageOption->GetFunction( 0 )->GetName(); +//modification{ + //ReadImage( inputImage, inputFile.c_str() ); + reader->SetFileName( inputFile.c_str() ); + inputImage = reader->GetOutput(); + inputImage->Update(); + inputImage->DisconnectPipeline(); +//}modification + } + else + { + if( verbose ) + { + std::cerr << "Input image not specified." << std::endl; + } + return EXIT_FAILURE; + } + + /** + * handle the mask image + */ + + bool isMaskImageSpecified = false; + + typename itk::ants::CommandLineParser::OptionType::Pointer maskImageOption = + parser->GetOption( "mask-image" ); + if( maskImageOption && maskImageOption->GetNumberOfFunctions() ) + { + std::string inputFile = maskImageOption->GetFunction( 0 )->GetName(); +//modification{ + //ReadImage( maskImage, inputFile.c_str() ); + typename ReaderType::Pointer maskreader = ReaderType::New(); + maskreader->SetFileName( inputFile.c_str() ); + maskImage = maskreader->GetOutput(); + maskImage->Update(); + maskImage->DisconnectPipeline(); +//}modification + + isMaskImageSpecified = true; + } + if( !maskImage ) + { + if( verbose ) + { + std::cout << "Mask not read. Using the entire image as the mask." << std::endl << std::endl; + } + maskImage = MaskImageType::New(); + maskImage->CopyInformation( inputImage ); + maskImage->SetRegions( inputImage->GetRequestedRegion() ); + maskImage->Allocate( false ); + maskImage->FillBuffer( itk::NumericTraits::OneValue() ); + } + + typename ImageType::Pointer weightImage = ITK_NULLPTR; + + typename itk::ants::CommandLineParser::OptionType::Pointer weightImageOption = + parser->GetOption( "weight-image" ); + if( weightImageOption && weightImageOption->GetNumberOfFunctions() ) + { + std::string inputFile = weightImageOption->GetFunction( 0 )->GetName(); +//modification{ + //ReadImage( weightImage, inputFile.c_str() ); + typename ReaderType::Pointer weightreader = ReaderType::New(); + weightreader->SetFileName( inputFile.c_str() ); + maskImage = weightreader->GetOutput(); + maskImage->Update(); + maskImage->DisconnectPipeline(); +//}modification + } + + /** + * convergence options + */ + typename itk::ants::CommandLineParser::OptionType::Pointer convergenceOption = + parser->GetOption( "convergence" ); + if( convergenceOption && convergenceOption->GetNumberOfFunctions() ) + { + if( convergenceOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + std::vector numIters = parser->ConvertVector( + convergenceOption->GetFunction( 0 )->GetParameter( 0 ) ); + typename CorrecterType::VariableSizeArrayType + maximumNumberOfIterations( numIters.size() ); + for( unsigned int d = 0; d < numIters.size(); d++ ) + { + maximumNumberOfIterations[d] = numIters[d]; + } + correcter->SetMaximumNumberOfIterations( maximumNumberOfIterations ); + + typename CorrecterType::ArrayType numberOfFittingLevels; + numberOfFittingLevels.Fill( numIters.size() ); + correcter->SetNumberOfFittingLevels( numberOfFittingLevels ); + correcter->SetConvergenceThreshold( 0.0 ); + } + if( convergenceOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + correcter->SetConvergenceThreshold( parser->Convert( + convergenceOption->GetFunction( 0 )->GetParameter( 1 ) ) ); + } + } + else // set default values + { + typename CorrecterType::VariableSizeArrayType + maximumNumberOfIterations( 4 ); + maximumNumberOfIterations.Fill( 50 ); + correcter->SetMaximumNumberOfIterations( maximumNumberOfIterations ); + correcter->SetNumberOfFittingLevels( 4 ); + correcter->SetConvergenceThreshold( 0.0 ); + } + + /** + * B-spline options -- we place this here to take care of the case where + * the user wants to specify things in terms of the spline distance. + */ + + typename ImageType::IndexType inputImageIndex = + inputImage->GetLargestPossibleRegion().GetIndex(); + typename ImageType::SizeType inputImageSize = + inputImage->GetLargestPossibleRegion().GetSize(); + + typename ImageType::PointType newOrigin = inputImage->GetOrigin(); + + typename itk::ants::CommandLineParser::OptionType::Pointer bsplineOption = + parser->GetOption( "bspline-fitting" ); + if( bsplineOption && bsplineOption->GetNumberOfFunctions() ) + { + if( bsplineOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + correcter->SetSplineOrder( parser->Convert( + bsplineOption->GetFunction( 0 )->GetParameter( 1 ) ) ); + } + if( bsplineOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + std::vector array = parser->ConvertVector( + bsplineOption->GetFunction( 0 )->GetParameter( 0 ) ); + typename CorrecterType::ArrayType numberOfControlPoints; + if( array.size() == 1 ) + { + // the user wants to specify things in terms of spline distance. + // 1. need to pad the images to get as close to possible to the + // requested domain size. + float splineDistance = array[0]; + + typename ImageType::SizeType originalImageSize = inputImage->GetLargestPossibleRegion().GetSize(); + + itk::Size lowerBound; + itk::Size upperBound; + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + float domain = static_cast( originalImageSize[d] - 1 ) * inputImage->GetSpacing()[d]; + unsigned int numberOfSpans = static_cast( + std::ceil( domain / splineDistance ) ); + unsigned long extraPadding = static_cast( ( numberOfSpans + * splineDistance + - domain ) / inputImage->GetSpacing()[d] + 0.5 ); + lowerBound[d] = static_cast( 0.5 * extraPadding ); + upperBound[d] = extraPadding - lowerBound[d]; + newOrigin[d] -= ( static_cast( lowerBound[d] ) + * inputImage->GetSpacing()[d] ); + numberOfControlPoints[d] = numberOfSpans + correcter->GetSplineOrder(); + } + + typedef itk::ConstantPadImageFilter PadderType; + typename PadderType::Pointer padder = PadderType::New(); + padder->SetInput( inputImage ); + padder->SetPadLowerBound( lowerBound ); + padder->SetPadUpperBound( upperBound ); + padder->SetConstant( 0 ); + padder->Update(); + + inputImage = padder->GetOutput(); + inputImage->DisconnectPipeline(); + + typedef itk::ConstantPadImageFilter MaskPadderType; + typename MaskPadderType::Pointer maskPadder = MaskPadderType::New(); + maskPadder->SetInput( maskImage ); + maskPadder->SetPadLowerBound( lowerBound ); + maskPadder->SetPadUpperBound( upperBound ); + maskPadder->SetConstant( 0 ); + maskPadder->Update(); + + maskImage = maskPadder->GetOutput(); + maskImage->DisconnectPipeline(); + + if( weightImage ) + { + typename PadderType::Pointer weightPadder = PadderType::New(); + weightPadder->SetInput( weightImage ); + weightPadder->SetPadLowerBound( lowerBound ); + weightPadder->SetPadUpperBound( upperBound ); + weightPadder->SetConstant( 0 ); + weightPadder->Update(); + + weightImage = weightPadder->GetOutput(); + weightImage->DisconnectPipeline(); + } + + if( verbose ) + { + std::cout << "Specified spline distance: " << splineDistance << "mm" << std::endl; + std::cout << " original image size: " << originalImageSize << std::endl; + std::cout << " padded image size: " << inputImage->GetLargestPossibleRegion().GetSize() << std::endl; + std::cout << " number of control points: " << numberOfControlPoints << std::endl; + std::cout << std::endl; + } + } + else if( array.size() == ImageDimension ) + { + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + numberOfControlPoints[d] = static_cast( array[d] ) + + correcter->GetSplineOrder(); + } + } + else + { + if( verbose ) + { + std::cerr << "Incorrect mesh resolution" << std::endl; + } + return EXIT_FAILURE; + } + correcter->SetNumberOfControlPoints( numberOfControlPoints ); + } + } + + typedef itk::ShrinkImageFilter ShrinkerType; + typename ShrinkerType::Pointer shrinker = ShrinkerType::New(); + shrinker->SetInput( inputImage ); + shrinker->SetShrinkFactors( 1 ); + + typedef itk::ShrinkImageFilter MaskShrinkerType; + typename MaskShrinkerType::Pointer maskshrinker = MaskShrinkerType::New(); + maskshrinker->SetInput( maskImage ); + maskshrinker->SetShrinkFactors( 1 ); + + typename itk::ants::CommandLineParser::OptionType::Pointer shrinkFactorOption = + parser->GetOption( "shrink-factor" ); + int shrinkFactor = 4; + if( shrinkFactorOption && shrinkFactorOption->GetNumberOfFunctions() ) + { + shrinkFactor = parser->Convert( shrinkFactorOption->GetFunction( 0 )->GetName() ); + } + shrinker->SetShrinkFactors( shrinkFactor ); + maskshrinker->SetShrinkFactors( shrinkFactor ); + if( ImageDimension == 4 ) + { + shrinker->SetShrinkFactor( 3, 1 ); + maskshrinker->SetShrinkFactor( 3, 1 ); + } + shrinker->Update(); + maskshrinker->Update(); + + itk::TimeProbe timer; + timer.Start(); + + correcter->SetInput( shrinker->GetOutput() ); + correcter->SetMaskImage( maskshrinker->GetOutput() ); + + typedef itk::ShrinkImageFilter WeightShrinkerType; + typename WeightShrinkerType::Pointer weightshrinker = WeightShrinkerType::New(); + if( weightImage ) + { + weightshrinker->SetInput( weightImage ); + weightshrinker->SetShrinkFactors( shrinker->GetShrinkFactors() ); + weightshrinker->Update(); + + correcter->SetConfidenceImage( weightshrinker->GetOutput() ); + } + + if( verbose ) + { + typedef CommandIterationUpdate CommandType; + typename CommandType::Pointer observer = CommandType::New(); + correcter->AddObserver( itk::IterationEvent(), observer ); + } + + /** + * histogram sharpening options + */ + typename itk::ants::CommandLineParser::OptionType::Pointer histOption = + parser->GetOption( "histogram-sharpening" ); + if( histOption && histOption->GetNumberOfFunctions() ) + { + if( histOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + correcter->SetBiasFieldFullWidthAtHalfMaximum( parser->Convert( + histOption->GetFunction( 0 )->GetParameter( 0 ) ) ); + } + if( histOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + correcter->SetWienerFilterNoise( parser->Convert( + histOption->GetFunction( 0 )->GetParameter( 1 ) ) ); + } + if( histOption->GetFunction( 0 )->GetNumberOfParameters() > 2 ) + { + correcter->SetNumberOfHistogramBins( parser->Convert( + histOption->GetFunction( 0 )->GetParameter( 2 ) ) ); + } + } + + try + { + // correcter->DebugOn(); + correcter->Update(); + } + catch( itk::ExceptionObject & e ) + { + if( verbose ) + { + std::cerr << "Exception caught: " << e << std::endl; + } + return EXIT_FAILURE; + } + + if( verbose ) + { + correcter->Print( std::cout, 3 ); + } + + timer.Stop(); + if( verbose ) + { + std::cout << "Elapsed time: " << timer.GetMean() << std::endl; + } + + /** + * output + */ + typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = + parser->GetOption( "output" ); + if( outputOption && outputOption->GetNumberOfFunctions() ) + { + /** + * Reconstruct the bias field at full image resolution. Divide + * the original input image by the bias field to get the final + * corrected image. + */ + typedef itk::BSplineControlPointImageFilter BSplinerType; + typename BSplinerType::Pointer bspliner = BSplinerType::New(); + bspliner->SetInput( correcter->GetLogBiasFieldControlPointLattice() ); + bspliner->SetSplineOrder( correcter->GetSplineOrder() ); + bspliner->SetSize( inputImage->GetLargestPossibleRegion().GetSize() ); + bspliner->SetOrigin( newOrigin ); + bspliner->SetDirection( inputImage->GetDirection() ); + bspliner->SetSpacing( inputImage->GetSpacing() ); + bspliner->Update(); + +//modification{ +// typename ImageType::Pointer logField = AllocImage( inputImage ); + typename ImageType::Pointer logField = ImageType::New(); + logField->SetOrigin( inputImage->GetOrigin() ); + logField->SetSpacing( inputImage->GetSpacing() ); + logField->SetRegions( inputImage->GetLargestPossibleRegion() ); + logField->SetDirection( inputImage->GetDirection() ); + logField->Allocate(); +//}modification + + itk::ImageRegionIterator ItB( + bspliner->GetOutput(), + bspliner->GetOutput()->GetLargestPossibleRegion() ); + itk::ImageRegionIterator ItF( logField, + logField->GetLargestPossibleRegion() ); + for( ItB.GoToBegin(), ItF.GoToBegin(); !ItB.IsAtEnd(); ++ItB, ++ItF ) + { + ItF.Set( ItB.Get()[0] ); + } + + typedef itk::ExpImageFilter ExpFilterType; + typename ExpFilterType::Pointer expFilter = ExpFilterType::New(); + expFilter->SetInput( logField ); + expFilter->Update(); + + typedef itk::DivideImageFilter DividerType; + typename DividerType::Pointer divider = DividerType::New(); + divider->SetInput1( inputImage ); + divider->SetInput2( expFilter->GetOutput() ); + divider->Update(); + + if( maskImageOption && maskImageOption->GetNumberOfFunctions() > 0 ) + { + itk::ImageRegionIteratorWithIndex ItD( divider->GetOutput(), + divider->GetOutput()->GetLargestPossibleRegion() ); + itk::ImageRegionIterator ItI( inputImage, + inputImage->GetLargestPossibleRegion() ); + for( ItD.GoToBegin(), ItI.GoToBegin(); !ItD.IsAtEnd(); ++ItD, ++ItI ) + { + if( maskImage->GetPixel( ItD.GetIndex() ) == itk::NumericTraits::ZeroValue() ) + { + ItD.Set( ItI.Get() ); + } + } + } + + bool doRescale = true; + + typename itk::ants::CommandLineParser::OptionType::Pointer rescaleOption = + parser->GetOption( "rescale-intensities" ); + if( ! isMaskImageSpecified || ( rescaleOption && rescaleOption->GetNumberOfFunctions() && + ! parser->Convert( rescaleOption->GetFunction()->GetName() ) ) ) + { + doRescale = false; + } + + if( doRescale ) + { + typedef itk::Image ShortImageType; + + typedef itk::BinaryThresholdImageFilter ThresholderType; + typename ThresholderType::Pointer thresholder = ThresholderType::New(); + thresholder->SetInsideValue( itk::NumericTraits::ZeroValue() ); + thresholder->SetOutsideValue( itk::NumericTraits::OneValue() ); + thresholder->SetLowerThreshold( itk::NumericTraits::ZeroValue() ); + thresholder->SetUpperThreshold( itk::NumericTraits::ZeroValue() ); + thresholder->SetInput( maskImage ); + + typedef itk::LabelStatisticsImageFilter StatsType; + typename StatsType::Pointer stats = StatsType::New(); + stats->SetInput( inputImage ); + stats->SetLabelInput( thresholder->GetOutput() ); + stats->UseHistogramsOff(); + stats->Update(); + + typedef typename StatsType::LabelPixelType StatsLabelType; + StatsLabelType maskLabel = itk::NumericTraits::OneValue(); + + RealType minOriginal = stats->GetMinimum( maskLabel ); + RealType maxOriginal = stats->GetMaximum( maskLabel ); + + typename StatsType::Pointer stats2 = StatsType::New(); + stats2->SetInput( divider->GetOutput() ); + stats2->SetLabelInput( thresholder->GetOutput() ); + stats2->UseHistogramsOff(); + stats2->Update(); + + RealType minBiasCorrected = stats2->GetMinimum( maskLabel ); + RealType maxBiasCorrected = stats2->GetMaximum( maskLabel ); + + RealType slope = ( maxOriginal - minOriginal ) / ( maxBiasCorrected - minBiasCorrected ); + + itk::ImageRegionIteratorWithIndex ItD( divider->GetOutput(), + divider->GetOutput()->GetLargestPossibleRegion() ); + for( ItD.GoToBegin(); !ItD.IsAtEnd(); ++ItD ) + { + if( maskImage->GetPixel( ItD.GetIndex() ) == maskLabel ) + { + RealType originalIntensity = ItD.Get(); + RealType rescaledIntensity = maxOriginal - slope * ( maxBiasCorrected - originalIntensity ); + ItD.Set( rescaledIntensity ); + } + } + } + + typename ImageType::RegionType inputRegion; + inputRegion.SetIndex( inputImageIndex ); + inputRegion.SetSize( inputImageSize ); + + typedef itk::ExtractImageFilter CropperType; + typename CropperType::Pointer cropper = CropperType::New(); + cropper->SetInput( divider->GetOutput() ); + cropper->SetExtractionRegion( inputRegion ); + cropper->SetDirectionCollapseToSubmatrix(); + cropper->Update(); + + typename CropperType::Pointer biasFieldCropper = CropperType::New(); + biasFieldCropper->SetInput( expFilter->GetOutput() ); + biasFieldCropper->SetExtractionRegion( inputRegion ); + biasFieldCropper->SetDirectionCollapseToSubmatrix(); + biasFieldCropper->Update(); + +//modification{ +/* + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() == 0 ) + { + WriteImage( cropper->GetOutput(), ( outputOption->GetFunction( 0 )->GetName() ).c_str() ); + } + else if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + WriteImage( cropper->GetOutput(), ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() ); + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + WriteImage( biasFieldCropper->GetOutput(), ( outputOption->GetFunction( 0 )->GetParameter( 1 ) ).c_str() ); + } + } + } +*/ + + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() == 0 ) + { + typedef itk::ImageFileWriter WriterType; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetInput( cropper->GetOutput() ); + writer->SetFileName( ( outputOption->GetFunction( 0 )->GetName() ).c_str() ); + writer->Update(); + } + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + typedef itk::ImageFileWriter WriterType; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetInput( cropper->GetOutput() ); + writer->SetFileName( ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() ); + writer->Update(); + } + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + typedef itk::ImageFileWriter WriterType; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetFileName( ( outputOption->GetFunction( 0 )->GetParameter( 1 ) ).c_str() ); + writer->SetInput( biasFieldCropper->GetOutput() ); + writer->Update(); + } + } +//}modification + + return EXIT_SUCCESS; +} + +void N4InitializeCommandLineOptions( itk::ants::CommandLineParser *parser ) +{ + typedef itk::ants::CommandLineParser::OptionType OptionType; + + { + std::string description = + std::string( "This option forces the image to be treated as a specified-" ) + + std::string( "dimensional image. If not specified, N4 tries to " ) + + std::string( "infer the dimensionality from the input image." ); + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "image-dimensionality" ); + option->SetShortName( 'd' ); + option->SetUsageOption( 0, "2/3/4" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "A scalar image is expected as input for bias correction. " ) + + std::string( "Since N4 log transforms the intensities, negative values " ) + + std::string( "or values close to zero should be processed prior to " ) + + std::string( "correction." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "input-image" ); + option->SetShortName( 'i' ); + option->SetUsageOption( 0, "inputImageFilename" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "If a mask image is specified, the final bias correction is " ) + + std::string( "only performed in the mask region. If a weight image is not " ) + + std::string( "specified, only intensity values inside the masked region are " ) + + std::string( "used during the execution of the algorithm. If a weight " ) + + std::string( "image is specified, only the non-zero weights are used in the " ) + + std::string( "execution of the algorithm although the mask region defines " ) + + std::string( "where bias correction is performed in the final output. " ) + + std::string( "Otherwise bias correction occurs over the entire image domain. " ) + + std::string( "See also the option description for the weight image. " ) + + std::string( "If a mask image is *not* specified then the entire image region " ) + + std::string( "will be used as the mask region. Note that this is different than " ) + + std::string( "the N3 implementation which uses the results of Otsu thresholding " ) + + std::string( "to define a mask. However, this leads to unknown anatomical regions being " ) + + std::string( "included and excluded during the bias correction." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "mask-image" ); + option->SetShortName( 'x' ); + option->SetUsageOption( 0, "maskImageFilename" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "At each iteration, a new intensity mapping is calculated " ) + + std::string( "and applied but there is nothing which constrains the " ) + + std::string( "new intensity range to be within certain values. The " ) + + std::string( "result is that the range can \"drift\" from the original " ) + + std::string( "at each iteration. This option rescales to the [min,max] " ) + + std::string( "range of the original image intensities within the user-specified mask." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "rescale-intensities" ); + option->SetShortName( 'r' ); + option->SetUsageOption( 0, "0/(1)" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "The weight image allows the user to perform a relative " ) + + std::string( "weighting of specific voxels during the B-spline fitting. " ) + + std::string( "For example, some studies have shown that N3 performed on " ) + + std::string( "white matter segmentations improves performance. If one " ) + + std::string( "has a spatial probability map of the white matter, one can " ) + + std::string( "use this map to weight the b-spline fitting towards those " ) + + std::string( "voxels which are more probabilistically classified as white " ) + + std::string( "matter. See also the option description for the mask image." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "weight-image" ); + option->SetUsageOption( 0, "weightImageFilename" ); + option->SetShortName( 'w' ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "Running N4 on large images can be time consuming. " ) + + std::string( "To lessen computation time, the input image can be resampled. " ) + + std::string( "The shrink factor, specified as a single integer, describes " ) + + std::string( "this resampling. Shrink factors <= 4 are commonly used." ) + + std::string( "Note that the shrink factor is only applied to the first two or " ) + + std::string( "three dimensions which we assume are spatial. " ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "shrink-factor" ); + option->SetShortName( 's' ); + option->SetUsageOption( 0, "1/2/3/(4)/..." ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "Convergence is determined by calculating the coefficient of " ) + + std::string( "variation between subsequent iterations. When this value " ) + + std::string( "is less than the specified threshold " ) + + std::string( "from the previous iteration or the maximum number of " ) + + std::string( "iterations is exceeded the program terminates. Multiple " ) + + std::string( "resolutions can be specified by using 'x' between the number " ) + + std::string( "of iterations at each resolution, e.g. 100x50x50." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "convergence" ); + option->SetShortName( 'c' ); + option->SetUsageOption( 0, "[,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "These options describe the b-spline fitting parameters. " ) + + std::string( "The initial b-spline mesh at the coarsest resolution is " ) + + std::string( "specified either as the number of elements in each dimension, " ) + + std::string( "e.g. 2x2x3 for 3-D images, or it can be specified as a " ) + + std::string( "single scalar parameter which describes the isotropic sizing " ) + + std::string( "of the mesh elements. The latter option is typically preferred. " ) + + std::string( "For each subsequent level, the spline distance decreases in " ) + + std::string( "half, or equivalently, the number of mesh elements doubles " ) + + std::string( "Cubic splines (order = 3) are typically used. The default setting " ) + + std::string( "is to employ a single mesh element over the entire domain, i.e., " ) + + std::string( "-b [1x1x1,3]." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "bspline-fitting" ); + option->SetShortName( 'b' ); + option->SetUsageOption( 0, "[splineDistance,]" ); + option->SetUsageOption( 1, "[initialMeshResolution,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "These options describe the histogram sharpening parameters, " ) + + std::string( "i.e. the deconvolution step parameters described in the " ) + + std::string( "original N3 algorithm. The default values have been shown " ) + + std::string( "to work fairly well." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "histogram-sharpening" ); + option->SetShortName( 't' ); + option->SetUsageOption( 0, "[,,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "The output consists of the bias corrected version of the " ) + + std::string( "input image. Optionally, one can also output the estimated " ) + + std::string( "bias field." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "output" ); + option->SetShortName( 'o' ); + option->SetUsageOption( 0, "correctedImage" ); + option->SetUsageOption( 1, "[correctedImage,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Get Version Information." ); + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "version" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Verbose output." ); + + OptionType::Pointer option = OptionType::New(); + option->SetShortName( 'v' ); + option->SetLongName( "verbose" ); + option->SetUsageOption( 0, "(0)/1" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Print the help menu (short version)." ); + + OptionType::Pointer option = OptionType::New(); + option->SetShortName( 'h' ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Print the help menu." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "help" ); + option->SetDescription( description ); + parser->AddOption( option ); + } +} + + +int main( int argc, char *argv[] ) +{ + +//modification{ +// entry point for the library; parameter 'args' is equivalent to 'argv' in (argc,argv) of commandline parameters to +// 'main()' +//int N4BiasFieldCorrection( std::vector args, std::ostream* /*out_stream = NULL */ ) +/*{ + // put the arguments coming in as 'args' into standard (argc,argv) format; + // 'args' doesn't have the command name as first, argument, so add it manually; + // 'args' may have adjacent arguments concatenated into one argument, + // which the parser should handle + args.insert( args.begin(), "N4BiasFieldCorrection" ); + + int argc = args.size(); + char* * argv = new char *[args.size() + 1]; + for( unsigned int i = 0; i < args.size(); ++i ) + { + // allocate space for the string plus a null character + argv[i] = new char[args[i].length() + 1]; + std::strncpy( argv[i], args[i].c_str(), args[i].length() ); + // place the null character in the end + argv[i][args[i].length()] = '\0'; + } + argv[argc] = ITK_NULLPTR; + // class to automatically cleanup argv upon destruction + class Cleanup_argv + { +public: + Cleanup_argv( char* * argv_, int argc_plus_one_ ) : argv( argv_ ), argc_plus_one( argc_plus_one_ ) + { + } + + ~Cleanup_argv() + { + for( unsigned int i = 0; i < argc_plus_one; ++i ) + { + delete[] argv[i]; + } + delete[] argv; + } + +private: + char* * argv; + unsigned int argc_plus_one; + }; + Cleanup_argv cleanup_argv( argv, argc + 1 ); + + // antscout->set_stream( out_stream ); +*/ +//}modification + + + itk::ants::CommandLineParser::Pointer parser = + itk::ants::CommandLineParser::New(); + + parser->SetCommand( argv[0] ); + + std::string commandDescription = + std::string( "N4 is a variant of the popular N3 (nonparameteric nonuniform " ) + + std::string( "normalization) retrospective bias correction algorithm. Based " ) + + std::string( "on the assumption that the corruption of the low frequency bias " ) + + std::string( "field can be modeled as a convolution of the intensity histogram " ) + + std::string( "by a Gaussian, the basic algorithmic protocol is to iterate " ) + + std::string( "between deconvolving the intensity histogram by a Gaussian, " ) + + std::string( "remapping the intensities, and then spatially smoothing this " ) + + std::string( "result by a B-spline modeling of the bias field itself. " ) + + std::string( "The modifications from and improvements obtained over " ) + + std::string( "the original N3 algorithm are described in the following paper: " ) + + std::string( "N. Tustison et al., N4ITK: Improved N3 Bias Correction, " ) + + std::string( "IEEE Transactions on Medical Imaging, 29(6):1310-1320, June 2010." ); + + parser->SetCommandDescription( commandDescription ); + N4InitializeCommandLineOptions( parser ); + + if( parser->Parse( argc, argv ) == EXIT_FAILURE ) + { + return EXIT_FAILURE; + } + + if( argc == 1 ) + { + parser->PrintMenu( std::cerr, 5, false ); + return EXIT_FAILURE; + } + else if( parser->GetOption( "help" )->GetFunction() && parser->Convert( parser->GetOption( "help" )->GetFunction()->GetName() ) ) + { + parser->PrintMenu( std::cout, 5, false ); + return EXIT_SUCCESS; + } + else if( parser->GetOption( 'h' )->GetFunction() && parser->Convert( parser->GetOption( 'h' )->GetFunction()->GetName() ) ) + { + parser->PrintMenu( std::cout, 5, true ); + return EXIT_SUCCESS; + } + +//modification{ +/* + // Show automatic version + itk::ants::CommandLineParser::OptionType::Pointer versionOption = parser->GetOption( "version" ); + if( versionOption && versionOption->GetNumberOfFunctions() ) + { + std::string versionFunction = versionOption->GetFunction( 0 )->GetName(); + ConvertToLowerCase( versionFunction ); + if( versionFunction.compare( "1" ) == 0 || versionFunction.compare( "true" ) == 0 ) + { + //Print Version Information + std::cout << ANTs::Version::ExtendedVersionString() << std::endl; + return EXIT_SUCCESS; + } + } +*/ +//}modification + + // Get dimensionality + unsigned int dimension = 3; + + itk::ants::CommandLineParser::OptionType::Pointer dimOption = + parser->GetOption( "image-dimensionality" ); + if( dimOption && dimOption->GetNumberOfFunctions() ) + { + dimension = parser->Convert( dimOption->GetFunction( 0 )->GetName() ); + } + else + { + // Read in the first intensity image to get the image dimension. + std::string filename; + + itk::ants::CommandLineParser::OptionType::Pointer imageOption = + parser->GetOption( "input-image" ); + if( imageOption && imageOption->GetNumberOfFunctions() > 0 ) + { + if( imageOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + filename = imageOption->GetFunction( 0 )->GetParameter( 0 ); + } + else + { + filename = imageOption->GetFunction( 0 )->GetName(); + } + } + else + { + std::cerr << "No input images were specified. Specify an input image" + << " with the -i option" << std::endl; + return EXIT_FAILURE; + } + itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO( + filename.c_str(), itk::ImageIOFactory::ReadMode ); + dimension = imageIO->GetNumberOfDimensions(); + } + + int returnValue = EXIT_FAILURE; + + switch( dimension ) + { + case 2: + { + returnValue = N4<2>( parser ); + } + break; + case 3: + { + returnValue = N4<3>( parser ); + } + break; + case 4: + { + returnValue = N4<4>( parser ); + } + break; + default: + std::cout << "Unsupported dimension" << std::endl; + return EXIT_FAILURE; + } + return returnValue; +} +//modification{ +//} // namespace ants +//}modification diff --git a/ThirdParty/ANTs/N4BiasFieldCorrection.cxx b/ThirdParty/ANTs/N4BiasFieldCorrection.cxx new file mode 100644 index 0000000..18e58df --- /dev/null +++ b/ThirdParty/ANTs/N4BiasFieldCorrection.cxx @@ -0,0 +1,947 @@ +#include "antsUtilities.h" +#include "antsAllocImage.h" +#include "antsCommandLineParser.h" + +#include "ReadWriteData.h" + +#include "itkBinaryThresholdImageFilter.h" +#include "itkBSplineControlPointImageFilter.h" +#include "itkConstantPadImageFilter.h" +#include "itkExpImageFilter.h" +#include "itkExtractImageFilter.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkLabelStatisticsImageFilter.h" +#include "itkN4BiasFieldCorrectionImageFilter.h" +#include "itkShrinkImageFilter.h" +#include "itkTimeProbe.h" + +#include +#include +#include + +#include "ANTsVersion.h" + +namespace ants +{ +template +class CommandIterationUpdate : public itk::Command +{ +public: + typedef CommandIterationUpdate Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + itkNewMacro( Self ); +protected: + CommandIterationUpdate() + { + }; +public: + + void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE + { + Execute( (const itk::Object *) caller, event); + } + + void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE + { + const TFilter * filter = + dynamic_cast( object ); + + if( typeid( event ) != typeid( itk::IterationEvent ) ) + { + return; + } + if( filter->GetElapsedIterations() == 1 ) + { + std::cout << "Current level = " << filter->GetCurrentLevel() + 1 + << std::endl; + } + std::cout << " Iteration " << filter->GetElapsedIterations() + << " (of " + << filter->GetMaximumNumberOfIterations()[filter->GetCurrentLevel()] + << "). "; + std::cout << " Current convergence value = " + << filter->GetCurrentConvergenceMeasurement() + << " (threshold = " << filter->GetConvergenceThreshold() + << ")" << std::endl; + } +}; + +template +int N4( itk::ants::CommandLineParser *parser ) +{ + typedef float RealType; + + typedef itk::Image ImageType; + typename ImageType::Pointer inputImage = ITK_NULLPTR; + + typedef itk::Image MaskImageType; + typename MaskImageType::Pointer maskImage = ITK_NULLPTR; + + bool verbose = false; + typename itk::ants::CommandLineParser::OptionType::Pointer verboseOption = + parser->GetOption( "verbose" ); + if( verboseOption && verboseOption->GetNumberOfFunctions() ) + { + verbose = parser->Convert( verboseOption->GetFunction( 0 )->GetName() ); + } + + if( verbose ) + { + std::cout << std::endl << "Running N4 for " + << ImageDimension << "-dimensional images." << std::endl << std::endl; + } + + typedef itk::N4BiasFieldCorrectionImageFilter CorrecterType; + typename CorrecterType::Pointer correcter = CorrecterType::New(); + typename itk::ants::CommandLineParser::OptionType::Pointer inputImageOption = + parser->GetOption( "input-image" ); + if( inputImageOption && inputImageOption->GetNumberOfFunctions() ) + { + std::string inputFile = inputImageOption->GetFunction( 0 )->GetName(); + ReadImage( inputImage, inputFile.c_str() ); + } + else + { + if( verbose ) + { + std::cerr << "Input image not specified." << std::endl; + } + return EXIT_FAILURE; + } + + /** + * handle the mask image + */ + + bool isMaskImageSpecified = false; + + typename itk::ants::CommandLineParser::OptionType::Pointer maskImageOption = + parser->GetOption( "mask-image" ); + if( maskImageOption && maskImageOption->GetNumberOfFunctions() ) + { + std::string inputFile = maskImageOption->GetFunction( 0 )->GetName(); + ReadImage( maskImage, inputFile.c_str() ); + + isMaskImageSpecified = true; + } + if( !maskImage ) + { + if( verbose ) + { + std::cout << "Mask not read. Using the entire image as the mask." << std::endl << std::endl; + } + maskImage = MaskImageType::New(); + maskImage->CopyInformation( inputImage ); + maskImage->SetRegions( inputImage->GetRequestedRegion() ); + maskImage->Allocate( false ); + maskImage->FillBuffer( itk::NumericTraits::OneValue() ); + } + + typename ImageType::Pointer weightImage = ITK_NULLPTR; + + typename itk::ants::CommandLineParser::OptionType::Pointer weightImageOption = + parser->GetOption( "weight-image" ); + if( weightImageOption && weightImageOption->GetNumberOfFunctions() ) + { + std::string inputFile = weightImageOption->GetFunction( 0 )->GetName(); + ReadImage( weightImage, inputFile.c_str() ); + } + + /** + * convergence options + */ + typename itk::ants::CommandLineParser::OptionType::Pointer convergenceOption = + parser->GetOption( "convergence" ); + if( convergenceOption && convergenceOption->GetNumberOfFunctions() ) + { + if( convergenceOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + std::vector numIters = parser->ConvertVector( + convergenceOption->GetFunction( 0 )->GetParameter( 0 ) ); + typename CorrecterType::VariableSizeArrayType + maximumNumberOfIterations( numIters.size() ); + for( unsigned int d = 0; d < numIters.size(); d++ ) + { + maximumNumberOfIterations[d] = numIters[d]; + } + correcter->SetMaximumNumberOfIterations( maximumNumberOfIterations ); + + typename CorrecterType::ArrayType numberOfFittingLevels; + numberOfFittingLevels.Fill( numIters.size() ); + correcter->SetNumberOfFittingLevels( numberOfFittingLevels ); + correcter->SetConvergenceThreshold( 0.0 ); + } + if( convergenceOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + correcter->SetConvergenceThreshold( parser->Convert( + convergenceOption->GetFunction( 0 )->GetParameter( 1 ) ) ); + } + } + else // set default values + { + typename CorrecterType::VariableSizeArrayType + maximumNumberOfIterations( 4 ); + maximumNumberOfIterations.Fill( 50 ); + correcter->SetMaximumNumberOfIterations( maximumNumberOfIterations ); + correcter->SetNumberOfFittingLevels( 4 ); + correcter->SetConvergenceThreshold( 0.0 ); + } + + /** + * B-spline options -- we place this here to take care of the case where + * the user wants to specify things in terms of the spline distance. + */ + + typename ImageType::IndexType inputImageIndex = + inputImage->GetLargestPossibleRegion().GetIndex(); + typename ImageType::SizeType inputImageSize = + inputImage->GetLargestPossibleRegion().GetSize(); + + typename ImageType::PointType newOrigin = inputImage->GetOrigin(); + + typename itk::ants::CommandLineParser::OptionType::Pointer bsplineOption = + parser->GetOption( "bspline-fitting" ); + if( bsplineOption && bsplineOption->GetNumberOfFunctions() ) + { + if( bsplineOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + correcter->SetSplineOrder( parser->Convert( + bsplineOption->GetFunction( 0 )->GetParameter( 1 ) ) ); + } + if( bsplineOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + std::vector array = parser->ConvertVector( + bsplineOption->GetFunction( 0 )->GetParameter( 0 ) ); + typename CorrecterType::ArrayType numberOfControlPoints; + if( array.size() == 1 ) + { + // the user wants to specify things in terms of spline distance. + // 1. need to pad the images to get as close to possible to the + // requested domain size. + float splineDistance = array[0]; + + typename ImageType::SizeType originalImageSize = inputImage->GetLargestPossibleRegion().GetSize(); + + itk::Size lowerBound; + itk::Size upperBound; + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + float domain = static_cast( originalImageSize[d] - 1 ) * inputImage->GetSpacing()[d]; + unsigned int numberOfSpans = static_cast( + std::ceil( domain / splineDistance ) ); + unsigned long extraPadding = static_cast( ( numberOfSpans + * splineDistance + - domain ) / inputImage->GetSpacing()[d] + 0.5 ); + lowerBound[d] = static_cast( 0.5 * extraPadding ); + upperBound[d] = extraPadding - lowerBound[d]; + newOrigin[d] -= ( static_cast( lowerBound[d] ) + * inputImage->GetSpacing()[d] ); + numberOfControlPoints[d] = numberOfSpans + correcter->GetSplineOrder(); + } + + typedef itk::ConstantPadImageFilter PadderType; + typename PadderType::Pointer padder = PadderType::New(); + padder->SetInput( inputImage ); + padder->SetPadLowerBound( lowerBound ); + padder->SetPadUpperBound( upperBound ); + padder->SetConstant( 0 ); + padder->Update(); + + inputImage = padder->GetOutput(); + inputImage->DisconnectPipeline(); + + typedef itk::ConstantPadImageFilter MaskPadderType; + typename MaskPadderType::Pointer maskPadder = MaskPadderType::New(); + maskPadder->SetInput( maskImage ); + maskPadder->SetPadLowerBound( lowerBound ); + maskPadder->SetPadUpperBound( upperBound ); + maskPadder->SetConstant( 0 ); + maskPadder->Update(); + + maskImage = maskPadder->GetOutput(); + maskImage->DisconnectPipeline(); + + if( weightImage ) + { + typename PadderType::Pointer weightPadder = PadderType::New(); + weightPadder->SetInput( weightImage ); + weightPadder->SetPadLowerBound( lowerBound ); + weightPadder->SetPadUpperBound( upperBound ); + weightPadder->SetConstant( 0 ); + weightPadder->Update(); + + weightImage = weightPadder->GetOutput(); + weightImage->DisconnectPipeline(); + } + + if( verbose ) + { + std::cout << "Specified spline distance: " << splineDistance << "mm" << std::endl; + std::cout << " original image size: " << originalImageSize << std::endl; + std::cout << " padded image size: " << inputImage->GetLargestPossibleRegion().GetSize() << std::endl; + std::cout << " number of control points: " << numberOfControlPoints << std::endl; + std::cout << std::endl; + } + } + else if( array.size() == ImageDimension ) + { + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + numberOfControlPoints[d] = static_cast( array[d] ) + + correcter->GetSplineOrder(); + } + } + else + { + if( verbose ) + { + std::cerr << "Incorrect mesh resolution" << std::endl; + } + return EXIT_FAILURE; + } + correcter->SetNumberOfControlPoints( numberOfControlPoints ); + } + } + + typedef itk::ShrinkImageFilter ShrinkerType; + typename ShrinkerType::Pointer shrinker = ShrinkerType::New(); + shrinker->SetInput( inputImage ); + shrinker->SetShrinkFactors( 1 ); + + typedef itk::ShrinkImageFilter MaskShrinkerType; + typename MaskShrinkerType::Pointer maskshrinker = MaskShrinkerType::New(); + maskshrinker->SetInput( maskImage ); + maskshrinker->SetShrinkFactors( 1 ); + + typename itk::ants::CommandLineParser::OptionType::Pointer shrinkFactorOption = + parser->GetOption( "shrink-factor" ); + int shrinkFactor = 4; + if( shrinkFactorOption && shrinkFactorOption->GetNumberOfFunctions() ) + { + shrinkFactor = parser->Convert( shrinkFactorOption->GetFunction( 0 )->GetName() ); + } + shrinker->SetShrinkFactors( shrinkFactor ); + maskshrinker->SetShrinkFactors( shrinkFactor ); + if( ImageDimension == 4 ) + { + shrinker->SetShrinkFactor( 3, 1 ); + maskshrinker->SetShrinkFactor( 3, 1 ); + } + shrinker->Update(); + maskshrinker->Update(); + + itk::TimeProbe timer; + timer.Start(); + + correcter->SetInput( shrinker->GetOutput() ); + correcter->SetMaskImage( maskshrinker->GetOutput() ); + + typedef itk::ShrinkImageFilter WeightShrinkerType; + typename WeightShrinkerType::Pointer weightshrinker = WeightShrinkerType::New(); + if( weightImage ) + { + weightshrinker->SetInput( weightImage ); + weightshrinker->SetShrinkFactors( shrinker->GetShrinkFactors() ); + weightshrinker->Update(); + + correcter->SetConfidenceImage( weightshrinker->GetOutput() ); + } + + if( verbose ) + { + typedef CommandIterationUpdate CommandType; + typename CommandType::Pointer observer = CommandType::New(); + correcter->AddObserver( itk::IterationEvent(), observer ); + } + + /** + * histogram sharpening options + */ + typename itk::ants::CommandLineParser::OptionType::Pointer histOption = + parser->GetOption( "histogram-sharpening" ); + if( histOption && histOption->GetNumberOfFunctions() ) + { + if( histOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + correcter->SetBiasFieldFullWidthAtHalfMaximum( parser->Convert( + histOption->GetFunction( 0 )->GetParameter( 0 ) ) ); + } + if( histOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + correcter->SetWienerFilterNoise( parser->Convert( + histOption->GetFunction( 0 )->GetParameter( 1 ) ) ); + } + if( histOption->GetFunction( 0 )->GetNumberOfParameters() > 2 ) + { + correcter->SetNumberOfHistogramBins( parser->Convert( + histOption->GetFunction( 0 )->GetParameter( 2 ) ) ); + } + } + + try + { + // correcter->DebugOn(); + correcter->Update(); + } + catch( itk::ExceptionObject & e ) + { + if( verbose ) + { + std::cerr << "Exception caught: " << e << std::endl; + } + return EXIT_FAILURE; + } + + if( verbose ) + { + correcter->Print( std::cout, 3 ); + } + + timer.Stop(); + if( verbose ) + { + std::cout << "Elapsed time: " << timer.GetMean() << std::endl; + } + + /** + * output + */ + typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = + parser->GetOption( "output" ); + if( outputOption && outputOption->GetNumberOfFunctions() ) + { + /** + * Reconstruct the bias field at full image resolution. Divide + * the original input image by the bias field to get the final + * corrected image. + */ + typedef itk::BSplineControlPointImageFilter BSplinerType; + typename BSplinerType::Pointer bspliner = BSplinerType::New(); + bspliner->SetInput( correcter->GetLogBiasFieldControlPointLattice() ); + bspliner->SetSplineOrder( correcter->GetSplineOrder() ); + bspliner->SetSize( inputImage->GetLargestPossibleRegion().GetSize() ); + bspliner->SetOrigin( newOrigin ); + bspliner->SetDirection( inputImage->GetDirection() ); + bspliner->SetSpacing( inputImage->GetSpacing() ); + bspliner->Update(); + + typename ImageType::Pointer logField = AllocImage( inputImage ); + + itk::ImageRegionIterator ItB( + bspliner->GetOutput(), + bspliner->GetOutput()->GetLargestPossibleRegion() ); + itk::ImageRegionIterator ItF( logField, + logField->GetLargestPossibleRegion() ); + for( ItB.GoToBegin(), ItF.GoToBegin(); !ItB.IsAtEnd(); ++ItB, ++ItF ) + { + ItF.Set( ItB.Get()[0] ); + } + + typedef itk::ExpImageFilter ExpFilterType; + typename ExpFilterType::Pointer expFilter = ExpFilterType::New(); + expFilter->SetInput( logField ); + expFilter->Update(); + + typedef itk::DivideImageFilter DividerType; + typename DividerType::Pointer divider = DividerType::New(); + divider->SetInput1( inputImage ); + divider->SetInput2( expFilter->GetOutput() ); + divider->Update(); + + if( maskImageOption && maskImageOption->GetNumberOfFunctions() > 0 ) + { + itk::ImageRegionIteratorWithIndex ItD( divider->GetOutput(), + divider->GetOutput()->GetLargestPossibleRegion() ); + itk::ImageRegionIterator ItI( inputImage, + inputImage->GetLargestPossibleRegion() ); + for( ItD.GoToBegin(), ItI.GoToBegin(); !ItD.IsAtEnd(); ++ItD, ++ItI ) + { + if( maskImage->GetPixel( ItD.GetIndex() ) == itk::NumericTraits::ZeroValue() ) + { + ItD.Set( ItI.Get() ); + } + } + } + + bool doRescale = true; + + typename itk::ants::CommandLineParser::OptionType::Pointer rescaleOption = + parser->GetOption( "rescale-intensities" ); + if( ! isMaskImageSpecified || ( rescaleOption && rescaleOption->GetNumberOfFunctions() && + ! parser->Convert( rescaleOption->GetFunction()->GetName() ) ) ) + { + doRescale = false; + } + + if( doRescale ) + { + typedef itk::Image ShortImageType; + + typedef itk::BinaryThresholdImageFilter ThresholderType; + typename ThresholderType::Pointer thresholder = ThresholderType::New(); + thresholder->SetInsideValue( itk::NumericTraits::ZeroValue() ); + thresholder->SetOutsideValue( itk::NumericTraits::OneValue() ); + thresholder->SetLowerThreshold( itk::NumericTraits::ZeroValue() ); + thresholder->SetUpperThreshold( itk::NumericTraits::ZeroValue() ); + thresholder->SetInput( maskImage ); + + typedef itk::LabelStatisticsImageFilter StatsType; + typename StatsType::Pointer stats = StatsType::New(); + stats->SetInput( inputImage ); + stats->SetLabelInput( thresholder->GetOutput() ); + stats->UseHistogramsOff(); + stats->Update(); + + typedef typename StatsType::LabelPixelType StatsLabelType; + StatsLabelType maskLabel = itk::NumericTraits::OneValue(); + + RealType minOriginal = stats->GetMinimum( maskLabel ); + RealType maxOriginal = stats->GetMaximum( maskLabel ); + + typename StatsType::Pointer stats2 = StatsType::New(); + stats2->SetInput( divider->GetOutput() ); + stats2->SetLabelInput( thresholder->GetOutput() ); + stats2->UseHistogramsOff(); + stats2->Update(); + + RealType minBiasCorrected = stats2->GetMinimum( maskLabel ); + RealType maxBiasCorrected = stats2->GetMaximum( maskLabel ); + + RealType slope = ( maxOriginal - minOriginal ) / ( maxBiasCorrected - minBiasCorrected ); + + itk::ImageRegionIteratorWithIndex ItD( divider->GetOutput(), + divider->GetOutput()->GetLargestPossibleRegion() ); + for( ItD.GoToBegin(); !ItD.IsAtEnd(); ++ItD ) + { + if( maskImage->GetPixel( ItD.GetIndex() ) == maskLabel ) + { + RealType originalIntensity = ItD.Get(); + RealType rescaledIntensity = maxOriginal - slope * ( maxBiasCorrected - originalIntensity ); + ItD.Set( rescaledIntensity ); + } + } + } + + typename ImageType::RegionType inputRegion; + inputRegion.SetIndex( inputImageIndex ); + inputRegion.SetSize( inputImageSize ); + + typedef itk::ExtractImageFilter CropperType; + typename CropperType::Pointer cropper = CropperType::New(); + cropper->SetInput( divider->GetOutput() ); + cropper->SetExtractionRegion( inputRegion ); + cropper->SetDirectionCollapseToSubmatrix(); + cropper->Update(); + + typename CropperType::Pointer biasFieldCropper = CropperType::New(); + biasFieldCropper->SetInput( expFilter->GetOutput() ); + biasFieldCropper->SetExtractionRegion( inputRegion ); + biasFieldCropper->SetDirectionCollapseToSubmatrix(); + biasFieldCropper->Update(); + + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() == 0 ) + { + WriteImage( cropper->GetOutput(), ( outputOption->GetFunction( 0 )->GetName() ).c_str() ); + } + else if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + WriteImage( cropper->GetOutput(), ( outputOption->GetFunction( 0 )->GetParameter( 0 ) ).c_str() ); + if( outputOption->GetFunction( 0 )->GetNumberOfParameters() > 1 ) + { + WriteImage( biasFieldCropper->GetOutput(), ( outputOption->GetFunction( 0 )->GetParameter( 1 ) ).c_str() ); + } + } + } + + return EXIT_SUCCESS; +} + +void N4InitializeCommandLineOptions( itk::ants::CommandLineParser *parser ) +{ + typedef itk::ants::CommandLineParser::OptionType OptionType; + + { + std::string description = + std::string( "This option forces the image to be treated as a specified-" ) + + std::string( "dimensional image. If not specified, N4 tries to " ) + + std::string( "infer the dimensionality from the input image." ); + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "image-dimensionality" ); + option->SetShortName( 'd' ); + option->SetUsageOption( 0, "2/3/4" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "A scalar image is expected as input for bias correction. " ) + + std::string( "Since N4 log transforms the intensities, negative values " ) + + std::string( "or values close to zero should be processed prior to " ) + + std::string( "correction." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "input-image" ); + option->SetShortName( 'i' ); + option->SetUsageOption( 0, "inputImageFilename" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "If a mask image is specified, the final bias correction is " ) + + std::string( "only performed in the mask region. If a weight image is not " ) + + std::string( "specified, only intensity values inside the masked region are " ) + + std::string( "used during the execution of the algorithm. If a weight " ) + + std::string( "image is specified, only the non-zero weights are used in the " ) + + std::string( "execution of the algorithm although the mask region defines " ) + + std::string( "where bias correction is performed in the final output. " ) + + std::string( "Otherwise bias correction occurs over the entire image domain. " ) + + std::string( "See also the option description for the weight image. " ) + + std::string( "If a mask image is *not* specified then the entire image region " ) + + std::string( "will be used as the mask region. Note that this is different than " ) + + std::string( "the N3 implementation which uses the results of Otsu thresholding " ) + + std::string( "to define a mask. However, this leads to unknown anatomical regions being " ) + + std::string( "included and excluded during the bias correction." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "mask-image" ); + option->SetShortName( 'x' ); + option->SetUsageOption( 0, "maskImageFilename" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "At each iteration, a new intensity mapping is calculated " ) + + std::string( "and applied but there is nothing which constrains the " ) + + std::string( "new intensity range to be within certain values. The " ) + + std::string( "result is that the range can \"drift\" from the original " ) + + std::string( "at each iteration. This option rescales to the [min,max] " ) + + std::string( "range of the original image intensities within the user-specified mask." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "rescale-intensities" ); + option->SetShortName( 'r' ); + option->SetUsageOption( 0, "0/(1)" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "The weight image allows the user to perform a relative " ) + + std::string( "weighting of specific voxels during the B-spline fitting. " ) + + std::string( "For example, some studies have shown that N3 performed on " ) + + std::string( "white matter segmentations improves performance. If one " ) + + std::string( "has a spatial probability map of the white matter, one can " ) + + std::string( "use this map to weight the b-spline fitting towards those " ) + + std::string( "voxels which are more probabilistically classified as white " ) + + std::string( "matter. See also the option description for the mask image." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "weight-image" ); + option->SetUsageOption( 0, "weightImageFilename" ); + option->SetShortName( 'w' ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "Running N4 on large images can be time consuming. " ) + + std::string( "To lessen computation time, the input image can be resampled. " ) + + std::string( "The shrink factor, specified as a single integer, describes " ) + + std::string( "this resampling. Shrink factors <= 4 are commonly used." ) + + std::string( "Note that the shrink factor is only applied to the first two or " ) + + std::string( "three dimensions which we assume are spatial. " ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "shrink-factor" ); + option->SetShortName( 's' ); + option->SetUsageOption( 0, "1/2/3/(4)/..." ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "Convergence is determined by calculating the coefficient of " ) + + std::string( "variation between subsequent iterations. When this value " ) + + std::string( "is less than the specified threshold " ) + + std::string( "from the previous iteration or the maximum number of " ) + + std::string( "iterations is exceeded the program terminates. Multiple " ) + + std::string( "resolutions can be specified by using 'x' between the number " ) + + std::string( "of iterations at each resolution, e.g. 100x50x50." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "convergence" ); + option->SetShortName( 'c' ); + option->SetUsageOption( 0, "[,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "These options describe the b-spline fitting parameters. " ) + + std::string( "The initial b-spline mesh at the coarsest resolution is " ) + + std::string( "specified either as the number of elements in each dimension, " ) + + std::string( "e.g. 2x2x3 for 3-D images, or it can be specified as a " ) + + std::string( "single scalar parameter which describes the isotropic sizing " ) + + std::string( "of the mesh elements. The latter option is typically preferred. " ) + + std::string( "For each subsequent level, the spline distance decreases in " ) + + std::string( "half, or equivalently, the number of mesh elements doubles " ) + + std::string( "Cubic splines (order = 3) are typically used. The default setting " ) + + std::string( "is to employ a single mesh element over the entire domain, i.e., " ) + + std::string( "-b [1x1x1,3]." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "bspline-fitting" ); + option->SetShortName( 'b' ); + option->SetUsageOption( 0, "[splineDistance,]" ); + option->SetUsageOption( 1, "[initialMeshResolution,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "These options describe the histogram sharpening parameters, " ) + + std::string( "i.e. the deconvolution step parameters described in the " ) + + std::string( "original N3 algorithm. The default values have been shown " ) + + std::string( "to work fairly well." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "histogram-sharpening" ); + option->SetShortName( 't' ); + option->SetUsageOption( 0, "[,,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = + std::string( "The output consists of the bias corrected version of the " ) + + std::string( "input image. Optionally, one can also output the estimated " ) + + std::string( "bias field." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "output" ); + option->SetShortName( 'o' ); + option->SetUsageOption( 0, "correctedImage" ); + option->SetUsageOption( 1, "[correctedImage,]" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Get Version Information." ); + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "version" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Verbose output." ); + + OptionType::Pointer option = OptionType::New(); + option->SetShortName( 'v' ); + option->SetLongName( "verbose" ); + option->SetUsageOption( 0, "(0)/1" ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Print the help menu (short version)." ); + + OptionType::Pointer option = OptionType::New(); + option->SetShortName( 'h' ); + option->SetDescription( description ); + parser->AddOption( option ); + } + + { + std::string description = std::string( "Print the help menu." ); + + OptionType::Pointer option = OptionType::New(); + option->SetLongName( "help" ); + option->SetDescription( description ); + parser->AddOption( option ); + } +} + +// entry point for the library; parameter 'args' is equivalent to 'argv' in (argc,argv) of commandline parameters to +// 'main()' +int N4BiasFieldCorrection( std::vector args, std::ostream* /*out_stream = NULL */ ) +{ + // put the arguments coming in as 'args' into standard (argc,argv) format; + // 'args' doesn't have the command name as first, argument, so add it manually; + // 'args' may have adjacent arguments concatenated into one argument, + // which the parser should handle + args.insert( args.begin(), "N4BiasFieldCorrection" ); + + int argc = args.size(); + char* * argv = new char *[args.size() + 1]; + for( unsigned int i = 0; i < args.size(); ++i ) + { + // allocate space for the string plus a null character + argv[i] = new char[args[i].length() + 1]; + std::strncpy( argv[i], args[i].c_str(), args[i].length() ); + // place the null character in the end + argv[i][args[i].length()] = '\0'; + } + argv[argc] = ITK_NULLPTR; + // class to automatically cleanup argv upon destruction + class Cleanup_argv + { +public: + Cleanup_argv( char* * argv_, int argc_plus_one_ ) : argv( argv_ ), argc_plus_one( argc_plus_one_ ) + { + } + + ~Cleanup_argv() + { + for( unsigned int i = 0; i < argc_plus_one; ++i ) + { + delete[] argv[i]; + } + delete[] argv; + } + +private: + char* * argv; + unsigned int argc_plus_one; + }; + Cleanup_argv cleanup_argv( argv, argc + 1 ); + + // antscout->set_stream( out_stream ); + + itk::ants::CommandLineParser::Pointer parser = + itk::ants::CommandLineParser::New(); + + parser->SetCommand( argv[0] ); + + std::string commandDescription = + std::string( "N4 is a variant of the popular N3 (nonparameteric nonuniform " ) + + std::string( "normalization) retrospective bias correction algorithm. Based " ) + + std::string( "on the assumption that the corruption of the low frequency bias " ) + + std::string( "field can be modeled as a convolution of the intensity histogram " ) + + std::string( "by a Gaussian, the basic algorithmic protocol is to iterate " ) + + std::string( "between deconvolving the intensity histogram by a Gaussian, " ) + + std::string( "remapping the intensities, and then spatially smoothing this " ) + + std::string( "result by a B-spline modeling of the bias field itself. " ) + + std::string( "The modifications from and improvements obtained over " ) + + std::string( "the original N3 algorithm are described in the following paper: " ) + + std::string( "N. Tustison et al., N4ITK: Improved N3 Bias Correction, " ) + + std::string( "IEEE Transactions on Medical Imaging, 29(6):1310-1320, June 2010." ); + + parser->SetCommandDescription( commandDescription ); + N4InitializeCommandLineOptions( parser ); + + if( parser->Parse( argc, argv ) == EXIT_FAILURE ) + { + return EXIT_FAILURE; + } + + if( argc == 1 ) + { + parser->PrintMenu( std::cerr, 5, false ); + return EXIT_FAILURE; + } + else if( parser->GetOption( "help" )->GetFunction() && parser->Convert( parser->GetOption( "help" )->GetFunction()->GetName() ) ) + { + parser->PrintMenu( std::cout, 5, false ); + return EXIT_SUCCESS; + } + else if( parser->GetOption( 'h' )->GetFunction() && parser->Convert( parser->GetOption( 'h' )->GetFunction()->GetName() ) ) + { + parser->PrintMenu( std::cout, 5, true ); + return EXIT_SUCCESS; + } + // Show automatic version + itk::ants::CommandLineParser::OptionType::Pointer versionOption = parser->GetOption( "version" ); + if( versionOption && versionOption->GetNumberOfFunctions() ) + { + std::string versionFunction = versionOption->GetFunction( 0 )->GetName(); + ConvertToLowerCase( versionFunction ); + if( versionFunction.compare( "1" ) == 0 || versionFunction.compare( "true" ) == 0 ) + { + //Print Version Information + std::cout << ANTs::Version::ExtendedVersionString() << std::endl; + return EXIT_SUCCESS; + } + } + // Get dimensionality + unsigned int dimension = 3; + + itk::ants::CommandLineParser::OptionType::Pointer dimOption = + parser->GetOption( "image-dimensionality" ); + if( dimOption && dimOption->GetNumberOfFunctions() ) + { + dimension = parser->Convert( dimOption->GetFunction( 0 )->GetName() ); + } + else + { + // Read in the first intensity image to get the image dimension. + std::string filename; + + itk::ants::CommandLineParser::OptionType::Pointer imageOption = + parser->GetOption( "input-image" ); + if( imageOption && imageOption->GetNumberOfFunctions() > 0 ) + { + if( imageOption->GetFunction( 0 )->GetNumberOfParameters() > 0 ) + { + filename = imageOption->GetFunction( 0 )->GetParameter( 0 ); + } + else + { + filename = imageOption->GetFunction( 0 )->GetName(); + } + } + else + { + std::cerr << "No input images were specified. Specify an input image" + << " with the -i option" << std::endl; + return EXIT_FAILURE; + } + itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO( + filename.c_str(), itk::ImageIOFactory::ReadMode ); + dimension = imageIO->GetNumberOfDimensions(); + } + + int returnValue = EXIT_FAILURE; + + switch( dimension ) + { + case 2: + { + returnValue = N4<2>( parser ); + } + break; + case 3: + { + returnValue = N4<3>( parser ); + } + break; + case 4: + { + returnValue = N4<4>( parser ); + } + break; + default: + std::cout << "Unsupported dimension" << std::endl; + return EXIT_FAILURE; + } + return returnValue; +} +} // namespace ants diff --git a/ThirdParty/ANTs/antsCommandLineOption.cxx b/ThirdParty/ANTs/antsCommandLineOption.cxx new file mode 100644 index 0000000..12b5c41 --- /dev/null +++ b/ThirdParty/ANTs/antsCommandLineOption.cxx @@ -0,0 +1,84 @@ +/*========================================================================= + + Program: Advanced Normalization Tools + + Copyright (c) ConsortiumOfANTS. All rights reserved. + See accompanying COPYING.txt or + https://github.com/stnava/ANTs/blob/master/ANTSCopyright.txt for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "antsCommandLineOption.h" + +namespace itk +{ +namespace ants +{ +CommandLineOption +::CommandLineOption() : m_ShortName( '\0' ), + m_LongName( "" ), + m_Description( "" ) +{ + this->m_OptionFunctions.clear(); + this->m_UsageOptions.clear(); +} + +void +CommandLineOption +::AddFunction( std::string functionString, char leftDelimiter, char rightDelimiter, unsigned int order ) +{ + OptionFunctionType::Pointer optionFunction = OptionFunctionType::New(); + + optionFunction->SetArgOrder( order ); + + std::string::size_type leftDelimiterPos = functionString.find( leftDelimiter ); + std::string::size_type rightDelimiterPos = functionString.find( rightDelimiter ); + + if( leftDelimiterPos == std::string::npos || + rightDelimiterPos == std::string::npos ) + { + optionFunction->SetName( functionString ); + this->m_OptionFunctions.push_front( optionFunction ); + } + else + { + OptionFunctionType::ParameterStackType parameters; + + optionFunction->SetName( functionString.substr( 0, leftDelimiterPos ) ); + + std::string::size_type leftPos = leftDelimiterPos; + std::string::size_type rightPos = functionString.find( ',', leftPos + 1 ); + while( rightPos != std::string::npos ) + { + parameters.push_back( functionString.substr( leftPos + 1, rightPos - leftPos - 1 ) ); + leftPos = rightPos; + rightPos = functionString.find( ',', leftPos + 1 ); + } + + rightPos = rightDelimiterPos; + parameters.push_back( functionString.substr( leftPos + 1, rightPos - leftPos - 1 ) ); + + optionFunction->SetParameters( parameters ); + + this->m_OptionFunctions.push_front( optionFunction ); + } + + this->Modified(); +} + +void +CommandLineOption +::SetUsageOption( unsigned int i, std::string usage ) +{ + if( i >= this->m_UsageOptions.size() ) + { + this->m_UsageOptions.resize( i + 1 ); + } + this->m_UsageOptions[i] = usage; + this->Modified(); +} +} // end namespace ants +} // end namespace itk diff --git a/ThirdParty/ANTs/antsCommandLineOption.h b/ThirdParty/ANTs/antsCommandLineOption.h new file mode 100644 index 0000000..f85a409 --- /dev/null +++ b/ThirdParty/ANTs/antsCommandLineOption.h @@ -0,0 +1,210 @@ +/*========================================================================= + + Program: Advanced Normalization Tools + + Copyright (c) ConsortiumOfANTS. All rights reserved. + See accompanying COPYING.txt or + https://github.com/stnava/ANTs/blob/master/ANTSCopyright.txt for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __antsCommandLineOption_h +#define __antsCommandLineOption_h + +#include "itkDataObject.h" +#include "itkObjectFactory.h" +#include "itkMacro.h" +#include "itkNumericTraits.h" + +#include +#include +#include + +namespace itk +{ +namespace ants +{ +/** \class CommandLineOption + \brief Simple data structure for holding command line options. + An option can have multiple values with each value holding 0 or more + parameters. E.g. suppose we were creating an image registration program + which has several transformation model options such as 'rigid', 'affine', + and 'deformable'. A instance of the command line option could have a + long name of "transformation", short name 't', and description + "Transformation model---rigid, affine, or deformable". The values for + this option would be "rigid", "affine", and "deformable". Each value + would then hold parameters that relate to that value. For example, a + possible subsection of the command line would be + + " --transformation rigid[parameter1,parameter2,etc.] + -m mutualinformation[parameter1] --optimization gradientdescent" +*/ + +class OptionFunction + : public DataObject +{ +public: + OptionFunction() : + m_Name( "" ), + m_ArgOrder( 0 ), + m_StageID( 0 ) + { + }; + ~OptionFunction() + { + }; + + typedef OptionFunction Self; + typedef DataObject Superclass; + typedef SmartPointer Pointer; + + itkNewMacro( Self ); + + itkTypeMacro( Option, DataObject ); + + typedef std::deque ParameterStackType; + + itkSetStringMacro( Name ); + itkGetStringMacro( Name ); + itkGetMacro( Name, std::string ); + + itkSetMacro( ArgOrder, unsigned int ); + itkGetConstMacro( ArgOrder, unsigned int ); + + itkSetMacro( StageID, unsigned int ); + itkGetConstMacro( StageID, unsigned int ); + + ParameterStackType GetParameters() + { + return this->m_Parameters; + } + + void SetParameters( ParameterStackType parameters ) + { + this->m_Parameters = parameters; + this->Modified(); + } + + std::string GetParameter( unsigned int i = 0 ) + { + if( i < this->m_Parameters.size() ) + { + return this->m_Parameters[i]; + } + else + { + std::string empty( "" ); + return empty; + } + } + + unsigned int GetNumberOfParameters() + { + return this->m_Parameters.size(); + } + +private: + std::string m_Name; + unsigned int m_ArgOrder; + unsigned int m_StageID; + ParameterStackType m_Parameters; +}; + +class CommandLineOption + : public DataObject +{ +public: + typedef CommandLineOption Self; + typedef DataObject Superclass; + typedef SmartPointer Pointer; + + itkNewMacro( Self ); + + itkTypeMacro( Option, DataObject ); + + typedef OptionFunction OptionFunctionType; + + typedef std::deque FunctionStackType; + typedef std::deque UsageOptionStackType; + + FunctionStackType GetFunctions() + { + return this->m_OptionFunctions; + } + + unsigned int GetNumberOfFunctions() + { + return this->m_OptionFunctions.size(); + } + + OptionFunction::Pointer GetFunction( unsigned int i = 0 ) + { + if( i < this->m_OptionFunctions.size() ) + { + return this->m_OptionFunctions[i]; + } + else + { + return ITK_NULLPTR; + } + } + + UsageOptionStackType GetUsageOptions() + { + return this->m_UsageOptions; + } + + unsigned int GetNumberOfUsageOptions() + { + return this->m_UsageOptions.size(); + } + + std::string GetUsageOption( unsigned int i = 0 ) + { + if( i < this->m_UsageOptions.size() ) + { + return this->m_UsageOptions[i]; + } + else + { + return std::string( "" ); + } + } + + itkSetMacro( ShortName, char ); + itkGetConstMacro( ShortName, char ); + + itkSetStringMacro( LongName ); + itkGetConstMacro( LongName, std::string ); + + itkSetStringMacro( Description ); + itkGetMacro( Description, std::string ); + + void AddFunction( std::string, char, char, unsigned int order = 0 ); + + void AddFunction( std::string s ) + { + this->AddFunction( s, '[', ']' ); + } + + void SetUsageOption( unsigned int, std::string ); + +protected: + CommandLineOption(); + virtual ~CommandLineOption() + { + }; +private: + char m_ShortName; + std::string m_LongName; + std::string m_Description; + UsageOptionStackType m_UsageOptions; + FunctionStackType m_OptionFunctions; +}; +} // end namespace ants +} // end namespace itk + +#endif diff --git a/ThirdParty/ANTs/antsCommandLineParser.cxx b/ThirdParty/ANTs/antsCommandLineParser.cxx new file mode 100644 index 0000000..8bb28b2 --- /dev/null +++ b/ThirdParty/ANTs/antsCommandLineParser.cxx @@ -0,0 +1,664 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "antsCommandLineParser.h" + +#include + +namespace itk +{ +namespace ants +{ + + std::string ConvertToHumanReadable(const std::string & input) + { + typedef std::map TypeMapType; + TypeMapType cnvtMap; + cnvtMap[typeid(signed char).name()] = "signed int"; + cnvtMap[typeid(unsigned char).name()] = "unsigned int"; + cnvtMap[typeid(signed short).name()] = "signed int"; + cnvtMap[typeid(unsigned short).name()] = "unsigned int"; + cnvtMap[typeid(signed int).name()] = "signed int"; + cnvtMap[typeid(unsigned int).name()] = "unsigned int"; + cnvtMap[typeid(signed long).name()] = "signed int"; + cnvtMap[typeid(unsigned long).name()] = "unsigned int"; + cnvtMap[typeid(float).name()] = "float"; + cnvtMap[typeid(double).name()] = "double"; + cnvtMap[typeid(std::string).name()] = "std::string"; + cnvtMap[typeid(char *).name()] = "char *"; + + TypeMapType::iterator mi=cnvtMap.find(input); + if ( mi == cnvtMap.end() ) + { + return std::string("Unmapped Type"); + } + return mi->second; + } +CommandLineParser +::CommandLineParser(): + m_LeftDelimiter ( '[' ), + m_RightDelimiter ( ']' ) +{ + this->m_Options.clear(); + this->m_Command.clear(); + this->m_CommandDescription.clear(); + this->m_UnknownOptions.clear(); +} + +void +CommandLineParser +::AddOption( OptionType::Pointer option ) +{ + if( ( option->GetShortName() != '\0' || + !this->GetOption( option->GetShortName() ) ) + || ( !option->GetLongName().empty() || + !this->GetOption( option->GetLongName() ) ) ) + { + this->m_Options.push_back( option ); + } + else + { + if( option->GetShortName() != '\0' && + this->GetOption( option->GetShortName() ) ) + { + itkWarningMacro( "Duplicate short option '-" + << option->GetShortName() << "'" ); + } + if( !( option->GetLongName().empty() ) && + this->GetOption( option->GetLongName() ) ) + { + itkWarningMacro( "Duplicate long option '--" + << option->GetLongName() << "'" ); + } + } +} + +bool +CommandLineParser +::starts_with(const std::string & s1, const std::string & s2) +{ + return s2.length() <= s1.length() && s1.compare(0, s2.length(), s2) == 0; +} + +int +CommandLineParser +::Parse( unsigned int argc, char * *argv ) +{ + std::vector arguments = + this->RegroupCommandLineArguments( argc, argv ); + + unsigned int n = 0; + unsigned int order = 0; + bool allFlagsAreValid = true; + + if ( arguments.size() > 1 ) + { + this->m_Command = arguments[n++]; + } + + while( n < arguments.size() ) + { + std::string argument = arguments[n++]; + std::string name; + + name.clear(); + if( starts_with( argument, "--" ) ) + { + name = argument.substr( 2, argument.length() - 1 ); + } + else if( starts_with( argument, "-" ) && !starts_with( argument, "--" ) ) + { + name = argument.substr( 1, 2 ); + } + if( atof( name.c_str() ) ) + { + continue; + } + if ( name.size() > 0 ) + { + allFlagsAreValid &= this->ValidateFlag( name ); + } + if( ( !( name.empty() ) ) && + ( name.size() > 0 ) + ) + { + OptionType::Pointer option = this->GetOption( name ); + if( !option ) + { + OptionType::Pointer unknownOption = OptionType::New(); + if( name.length() > 1 ) + { + unknownOption->SetLongName( name ); + } + else + { + unknownOption->SetShortName( name.at( 0 ) ); + } + if( n == arguments.size() ) + { + unknownOption->AddFunction( "1", + this->m_LeftDelimiter, this->m_RightDelimiter, order++ ); + } + else + { + for( unsigned int m = n; m < arguments.size(); m++ ) + { + std::string function = arguments[m]; + if( !starts_with( function, "-" ) ) + { + unknownOption->AddFunction( function, + this->m_LeftDelimiter, this->m_RightDelimiter, order++ ); + } + else + { + if( m == n ) + { + unknownOption->AddFunction( "1", + this->m_LeftDelimiter, this->m_RightDelimiter, order++ ); + } + n = m; + break; + } + } + } + this->m_UnknownOptions.push_back( unknownOption ); + } + else // the option exists + { + if( n == arguments.size() ) + { + option->AddFunction( "1", + this->m_LeftDelimiter, this->m_RightDelimiter, order++ ); + } + else + { + for( unsigned int m = n; m < arguments.size(); m++ ) + { + std::string function = arguments[m]; + if( !starts_with( function, "-" ) || atof( function.c_str() ) ) + { + option->AddFunction( function, + this->m_LeftDelimiter, this->m_RightDelimiter, order++ ); + } + else + { + if( m == n ) + { + option->AddFunction( "1", + this->m_LeftDelimiter, this->m_RightDelimiter, order++ ); + } + n = m; + break; + } + } + } + } + } + } + if( ! allFlagsAreValid ) + { + std::cerr << "ERROR: Invalid command line flags found! Aborting execution." << std::endl; + return EXIT_FAILURE; + } + this->AssignStages(); + return EXIT_SUCCESS; +} + +std::vector +CommandLineParser +::RegroupCommandLineArguments( unsigned int argc, char * *argv ) +{ + /** + * Inclusion of this function allows the user to use spaces inside + * the left and right delimiters. Also replace other left and right + * delimiters. + */ + std::vector arguments; + + std::string currentArg( "" ); + bool isArgOpen = false; + + for( unsigned int n = 0; n < argc; n++ ) + { + std::string a( argv[n] ); + + if( n == 0 ) + { + arguments.push_back( a ); + } + else + { + // replace left delimiters + std::replace( a.begin(), a.begin()+1, '{', '[' ); + std::replace( a.begin(), a.begin()+1, '(', '[' ); + std::replace( a.begin(), a.begin()+1, '<', '[' ); + + // replace right delimiters + std::replace( a.end()-1, a.end(), '}', ']' ); + std::replace( a.end()-1, a.end(), ')', ']' ); + std::replace( a.end()-1, a.end(), '>', ']' ); + + if( isArgOpen ) + { + std::size_t leftDelimiterPosition = a.find( this->m_LeftDelimiter ); + if( leftDelimiterPosition != std::string::npos ) + { + itkExceptionMacro( "Incorrect command line specification. Missing leftDelimiterPosition? " << a ); + } + + std::size_t rightDelimiterPosition = a.find( this->m_RightDelimiter ); + if( rightDelimiterPosition != std::string::npos ) + { + if( rightDelimiterPosition < a.length() - 1 ) + { + itkExceptionMacro( "Incorrect command line specification. Missing rightDelimiterPosition? " << a ); + } + else + { + currentArg += a; + arguments.push_back( currentArg ); + currentArg.clear(); + isArgOpen = false; + } + } + else + { + currentArg += a; + } + } + else + { + std::size_t leftDelimiterPosition = a.find( this->m_LeftDelimiter ); + std::size_t rightDelimiterPosition = a.find( this->m_RightDelimiter ); + + if( leftDelimiterPosition == std::string::npos ) + { + if( rightDelimiterPosition == std::string::npos ) + { + currentArg += a; + arguments.push_back( currentArg ); + currentArg.clear(); + } + else + { + itkExceptionMacro( "Incorrect command line specification. " << a); + } + } + else if( leftDelimiterPosition != std::string::npos && + rightDelimiterPosition != std::string::npos && + leftDelimiterPosition < rightDelimiterPosition ) + { + if( rightDelimiterPosition < a.length() - 1 ) + { + itkExceptionMacro( "Incorrect command line specification. " << a ); + } + currentArg += a; + arguments.push_back( currentArg ); + currentArg.clear(); + isArgOpen = false; + } + else if( rightDelimiterPosition == std::string::npos && + leftDelimiterPosition != std::string::npos ) + { + currentArg += a; + isArgOpen = true; + } + } + } + } + + return arguments; +} + +CommandLineParser::OptionType::Pointer +CommandLineParser +::GetOption( std::string name ) +{ + if( name.length() == 1 ) + { + return this->GetOption( name.at( 0 ) ); + } + + OptionListType::iterator it; + for( it = this->m_Options.begin(); it != this->m_Options.end(); ++it ) + { + if( name.compare( (*it)->GetLongName() ) == 0 ) + { + return *it; + } + } + return ITK_NULLPTR; +} + +CommandLineParser::OptionType::Pointer +CommandLineParser +::GetOption( char name ) +{ + OptionListType::iterator it; + + for( it = this->m_Options.begin(); it != this->m_Options.end(); ++it ) + { + if( name == (*it)->GetShortName() ) + { + return *it; + } + } + return ITK_NULLPTR; +} + +bool +CommandLineParser:: +ValidateFlag(const std::string & currentFlag) +{ + bool validFlagFound = false; + for( OptionListType::const_iterator it = this->m_Options.begin(); it != this->m_Options.end(); ++it ) + { + const char shortName = (*it)->GetShortName(); + const std::string longName = (*it)->GetLongName(); + if( ( ( currentFlag.size() == 1 ) && + ( shortName == currentFlag[0] ) ) + || ( longName == currentFlag ) ) + { + validFlagFound = true; + } + } + + if ( ( ! validFlagFound ) && ( currentFlag.size() > 0 ) && ( !atof( currentFlag.c_str() ) ) ) + { + std::cout << "ERROR: Invalid flag provided " << currentFlag << std::endl; + } + return validFlagFound; +} + +void +CommandLineParser +::PrintMenu( std::ostream& os, Indent indent, bool printShortVersion ) const +{ + os << std::endl; + os << "COMMAND: " << std::endl; + os << indent << this->m_Command << std::endl; + if( !this->m_CommandDescription.empty() && !printShortVersion ) + { + std::stringstream ss1; + ss1 << indent << indent; + + std::stringstream ss2; + ss2 << this->m_CommandDescription; + + std::string description = this->BreakUpStringIntoNewLines( + ss2.str(), ss1.str(), 80 ); + + os << indent << indent << description << std::endl; + } + os << std::endl; + os << "OPTIONS: " << std::endl; + + OptionListType::const_iterator it; + for( it = this->m_Options.begin(); it != this->m_Options.end(); ++it ) + { + os << indent; + std::stringstream ss; + ss << indent; + + if( (*it)->GetShortName() != '\0' ) + { + os << "-" << (*it)->GetShortName(); + ss << Indent( 2 ); + if( !( (*it)->GetLongName() ).empty() ) + { + os << ", " << "--" << (*it)->GetLongName() << " " << std::flush; + ss << Indent( 5 + ( (*it)->GetLongName() ).length() ); + } + else + { + os << " " << std::flush; + ss << Indent( 1 ); + } + } + else + { + os << "--" << (*it)->GetLongName() << " " << std::flush; + ss << Indent( 3 + ( (*it)->GetLongName() ).length() ); + } + if( (*it)->GetNumberOfUsageOptions() > 0 ) + { + os << (*it)->GetUsageOption( 0 ) << std::endl; + for( unsigned int i = 1; i < (*it)->GetNumberOfUsageOptions(); i++ ) + { + os << ss.str() << (*it)->GetUsageOption( i ) << std::endl; + } + } + else + { + os << std::endl; + } + + if( !( (*it)->GetDescription().empty() ) && !printShortVersion ) + { + std::stringstream ss1; + ss1 << indent << indent; + + std::stringstream ss2; + ss2 << (*it)->GetDescription(); + + std::string description = this->BreakUpStringIntoNewLines( ss2.str(), ss1.str(), 80 ); + + os << indent << indent << description << std::endl; + } + if( !printShortVersion ) + { + if( (*it)->GetFunctions().size() == 1 ) + { + os << indent << indent << ": " << (*it)->GetFunction( 0 )->GetName(); + if( (*it)->GetFunction( 0 )->GetParameters().size() > 0 ) + { + os << "["; + if( (*it)->GetFunction( 0 )->GetParameters().size() == 1 ) + { + os << (*it)->GetFunction( 0 )->GetParameter( 0 ); + } + else + { + for( unsigned int i = 0; i < (*it)->GetFunction( 0 )->GetParameters().size() - 1; i++ ) + { + os << (*it)->GetFunction( 0 )->GetParameter( i ) << ","; + } + os << (*it)->GetFunction( 0 )->GetParameter( (*it)->GetFunction( 0 )->GetParameters().size() - 1 ); + } + os << "]"; + } + os << std::endl; + } + else if( (*it)->GetFunctions().size() > 1 ) + { + os << indent << indent << ": "; + for( unsigned int n = 0; n < (*it)->GetFunctions().size() - 1; n++ ) + { + os << (*it)->GetFunction( n )->GetName(); + if( (*it)->GetFunction( n )->GetParameters().size() > 0 ) + { + os << "["; + if( (*it)->GetFunction( n )->GetParameters().size() == 1 ) + { + os << (*it)->GetFunction( n )->GetParameter( 0 ) << "], "; + } + else + { + for( unsigned int i = 0; i < (*it)->GetFunction( n )->GetParameters().size() - 1; i++ ) + { + os << (*it)->GetFunction( n )->GetParameter( i ) << ","; + } + os + << (*it)->GetFunction( n )->GetParameter( (*it)->GetFunction( n )->GetParameters().size() + - 1 ) << "], "; + } + } + else + { + os << ", "; + } + } + + unsigned int nn = (*it)->GetFunctions().size() - 1; + + os << (*it)->GetFunction( nn )->GetName(); + if( (*it)->GetFunction( nn )->GetParameters().size() > 0 ) + { + os << "["; + if( (*it)->GetFunction( nn )->GetParameters().size() == 1 ) + { + os << (*it)->GetFunction( nn )->GetParameter( 0 ) << "]"; + } + else + { + for( unsigned int i = 0; i < (*it)->GetFunction( nn )->GetParameters().size() - 1; i++ ) + { + os << (*it)->GetFunction( nn )->GetParameter( i ) << ","; + } + os << (*it)->GetFunction( nn )->GetParameter( (*it)->GetFunction( nn )->GetParameters().size() - 1 ) << "]"; + } + } + } + os << std::endl; + } + } +} + +std::string +CommandLineParser +::BreakUpStringIntoNewLines( std::string longString, + std::string indentString, unsigned int numberOfCharactersPerLine ) const +{ + std::vector tokens; + + this->TokenizeString( longString, tokens, " " ); + + std::string newString( "" ); + unsigned int currentTokenId = 0; + unsigned int currentLineLength = 0; + while( currentTokenId < tokens.size() ) + { + if( tokens[currentTokenId].length() >= numberOfCharactersPerLine ) + { + newString += ( std::string( "\n" ) + tokens[currentTokenId] + + std::string( "\n" ) ); + currentTokenId++; + currentLineLength = 0; + } + else if( currentTokenId < tokens.size() && currentLineLength + + tokens[currentTokenId].length() > numberOfCharactersPerLine ) + { + newString += ( std::string( "\n" ) + indentString ); + currentLineLength = 0; + } + else + { + newString += ( tokens[currentTokenId] + std::string( " " ) ); + currentLineLength += ( tokens[currentTokenId].length() + 1 ); + currentTokenId++; + } + } + + return newString; +} + +void +CommandLineParser +::TokenizeString( std::string str, std::vector & tokens, + std::string delimiters ) const +{ + // Skip delimiters at beginning. + std::string::size_type lastPos = str.find_first_not_of( delimiters, 0 ); + // Find first "non-delimiter". + std::string::size_type pos = str.find_first_of( delimiters, lastPos ); + + while( std::string::npos != pos || std::string::npos != lastPos ) + { + // Found a token, add it to the vector. + tokens.push_back( str.substr( lastPos, pos - lastPos ) ); + // Skip delimiters. Note the "not_of" + lastPos = str.find_first_not_of( delimiters, pos ); + // Find next "non-delimiter" + pos = str.find_first_of( delimiters, lastPos ); + } +} + +void +CommandLineParser +::AssignStages() +{ + OptionListType::const_iterator it; + + for( it = this->m_Options.begin(); it != this->m_Options.end(); ++it ) + { + typedef OptionType::FunctionStackType OptionFunctionStackType; + OptionFunctionStackType functions = (*it)->GetFunctions(); + + OptionFunctionStackType::const_iterator it2; + + unsigned int previousOrder = 0; + unsigned int currentOrder = 0; + for( it2 = functions.begin(); it2 != functions.end(); ++it2 ) + { + if( it2 == functions.begin() ) + { + previousOrder = (*it2)->GetArgOrder(); + (*it2)->SetStageID( 0 ); + } + else + { + currentOrder = (*it2)->GetArgOrder(); + if( previousOrder == currentOrder + 1 ) + { + (*it2)->SetStageID( functions[it2 - functions.begin() - 1]->GetStageID() ); + } + else + { + (*it2)->SetStageID( functions[it2 - functions.begin() - 1]->GetStageID() + 1 ); + } + previousOrder = currentOrder; + } + } + } +} + +/** + * Standard "PrintSelf" method + */ +void +CommandLineParser +::PrintSelf( std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf( os, indent ); + + os << indent << "Command: " << this->m_Command << std::endl; + os << indent << "Options: " << std::endl; + + OptionListType::const_iterator it; + for( it = this->m_Options.begin(); it != this->m_Options.end(); ++it ) + { + (*it)->Print( os, indent ); + } + + if( this->m_UnknownOptions.size() ) + { + os << indent << "Unknown Options: " << std::endl; + OptionListType::const_iterator its; + for( its = this->m_UnknownOptions.begin(); + its != this->m_UnknownOptions.end(); ++its ) + { + (*its)->Print( os, indent ); + } + } +} +} // end namespace ants +} // end namespace itk diff --git a/ThirdParty/ANTs/antsCommandLineParser.h b/ThirdParty/ANTs/antsCommandLineParser.h new file mode 100644 index 0000000..1fd3133 --- /dev/null +++ b/ThirdParty/ANTs/antsCommandLineParser.h @@ -0,0 +1,184 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __antsCommandLineParser_h +#define __antsCommandLineParser_h + +#include "antsCommandLineOption.h" + +#include "itkDataObject.h" +#include "itkObjectFactory.h" +#include "itkMacro.h" +#include "itkNumericTraits.h" + +#include +#include +#include +#include +#include + +#include + +namespace itk +{ +namespace ants +{ + /** + * A untilty function to convert internal typeid(???).name() to + * the human readable equivalent format. + */ + extern std::string ConvertToHumanReadable(const std::string & input); + +/** \class CommandLineParser + \brief Simple command line parser. + \par + Parses the standard ( argc, argv ) variables which are stored + as options in the helper class antsCommandLineOption. Also contains + routines for converting types including std::vectors using 'x' as a + delimiter. For example, I can specify the 3-element std::vector + {10, 20, 30} as "10x20x30". +*/ + +class CommandLineParser + : public DataObject +{ +public: + /** Standard class typedefs. */ + typedef CommandLineParser Self; + typedef DataObject Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro( Self ); + + /** Run-time type information (and related methods). */ + itkTypeMacro( CommandLineParser, DataObject ); + + typedef CommandLineOption OptionType; + typedef std::list OptionListType; + typedef std::list StringListType; + + /** + * Interface routines + */ + + OptionType::Pointer GetOption( char ); + + OptionType::Pointer GetOption( std::string ); + + bool starts_with( const std::string &, const std::string & ); + int Parse( unsigned int, char * * ); + bool ValidateFlag(const std::string & currentFlag); + + void AddOption( OptionType::Pointer ); + + void PrintMenu( std::ostream& os, Indent indent, bool printShortVersion = false ) const; + + itkSetStringMacro( Command ); + itkGetStringMacro( Command ); + + itkSetStringMacro( CommandDescription ); + itkGetStringMacro( CommandDescription ); + + OptionListType GetOptions() const + { + return this->m_Options; + } + + OptionListType GetUnknownOptions() const + { + return this->m_UnknownOptions; + } + + /** + * This feature is designed for a more advanced command line usage + * where multiple option values are used per stage (e.g. + * antsRegistration). Multiple option value are considered to be of + * the same stage if they are situated adjacently on the command line. + */ + void AssignStages(); + + + template + TValue Convert( std::string optionString ) const + { + //Strip whitespace at end + optionString.erase(optionString.find_last_not_of(" \n\r\t")+1); + TValue value; + std::istringstream iss( optionString ); + if (!(iss >> value) //Conversion did not fail + || !( iss.peek() == EOF ) // All content parsed + ) + { + std::string internalTypeName( typeid(value).name() ); + itkExceptionMacro( "ERROR: Parse error occured during command line argument processing\n" + << "ERROR: Unable to convert '" << optionString + << "' to type '" << internalTypeName << "' as " << ConvertToHumanReadable(internalTypeName) << std::endl); + } + return value; + } + + template + std::vector ConvertVector( std::string optionString ) const + { + //Strip whitespace at end + optionString.erase(optionString.find_last_not_of(" \n\r\t")+1); + + std::vector optionElementString; + std::istringstream f(optionString); + std::string s; + while( std::getline(f, s, 'x')) + { + optionElementString.push_back(s); + } + + std::vector< TValue > values; + for ( std::vector< std::string >::const_iterator oESit = optionElementString.begin(); + oESit != optionElementString.end(); ++oESit) + { + const TValue & value = this->Convert( *oESit ); + values.push_back ( value ); + } + return values; + } + +protected: + CommandLineParser(); + virtual ~CommandLineParser() + { + } + + void PrintSelf( std::ostream& os, Indent indent ) const ITK_OVERRIDE; + +private: + CommandLineParser( const Self & ); // purposely not implemented + void operator=( const Self & ); // purposely not implemented + + std::vector RegroupCommandLineArguments( unsigned int, char * * ); + + std::string BreakUpStringIntoNewLines( std::string, const std::string, unsigned int ) const; + + void TokenizeString( std::string, std::vector &, std::string ) const; + + OptionListType m_Options; + std::string m_Command; + std::string m_CommandDescription; + OptionListType m_UnknownOptions; + + char m_LeftDelimiter; + char m_RightDelimiter; +}; +} // end namespace ants +} // end namespace itk + +#endif diff --git a/ThirdParty/ITK/LICENSE b/ThirdParty/ITK/LICENSE deleted file mode 100644 index d645695..0000000 --- a/ThirdParty/ITK/LICENSE +++ /dev/null @@ -1,202 +0,0 @@ - - Apache License - Version 2.0, January 2004 - http://www.apache.org/licenses/ - - TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION - - 1. Definitions. - - "License" shall mean the terms and conditions for use, reproduction, - and distribution as defined by Sections 1 through 9 of this document. - - "Licensor" shall mean the copyright owner or entity authorized by - the copyright owner that is granting the License. - - "Legal Entity" shall mean the union of the acting entity and all - other entities that control, are controlled by, or are under common - control with that entity. For the purposes of this definition, - "control" means (i) the power, direct or indirect, to cause the - direction or management of such entity, whether by contract or - otherwise, or (ii) ownership of fifty percent (50%) or more of the - outstanding shares, or (iii) beneficial ownership of such entity. - - "You" (or "Your") shall mean an individual or Legal Entity - exercising permissions granted by this License. - - "Source" form shall mean the preferred form for making modifications, - including but not limited to software source code, documentation - source, and configuration files. - - "Object" form shall mean any form resulting from mechanical - transformation or translation of a Source form, including but - not limited to compiled object code, generated documentation, - and conversions to other media types. - - "Work" shall mean the work of authorship, whether in Source or - Object form, made available under the License, as indicated by a - copyright notice that is included in or attached to the work - (an example is provided in the Appendix below). - - "Derivative Works" shall mean any work, whether in Source or Object - form, that is based on (or derived from) the Work and for which the - editorial revisions, annotations, elaborations, or other modifications - represent, as a whole, an original work of authorship. For the purposes - of this License, Derivative Works shall not include works that remain - separable from, or merely link (or bind by name) to the interfaces of, - the Work and Derivative Works thereof. - - "Contribution" shall mean any work of authorship, including - the original version of the Work and any modifications or additions - to that Work or Derivative Works thereof, that is intentionally - submitted to Licensor for inclusion in the Work by the copyright owner - or by an individual or Legal Entity authorized to submit on behalf of - the copyright owner. For the purposes of this definition, "submitted" - means any form of electronic, verbal, or written communication sent - to the Licensor or its representatives, including but not limited to - communication on electronic mailing lists, source code control systems, - and issue tracking systems that are managed by, or on behalf of, the - Licensor for the purpose of discussing and improving the Work, but - excluding communication that is conspicuously marked or otherwise - designated in writing by the copyright owner as "Not a Contribution." - - "Contributor" shall mean Licensor and any individual or Legal Entity - on behalf of whom a Contribution has been received by Licensor and - subsequently incorporated within the Work. - - 2. Grant of Copyright License. Subject to the terms and conditions of - this License, each Contributor hereby grants to You a perpetual, - worldwide, non-exclusive, no-charge, royalty-free, irrevocable - copyright license to reproduce, prepare Derivative Works of, - publicly display, publicly perform, sublicense, and distribute the - Work and such Derivative Works in Source or Object form. - - 3. Grant of Patent License. Subject to the terms and conditions of - this License, each Contributor hereby grants to You a perpetual, - worldwide, non-exclusive, no-charge, royalty-free, irrevocable - (except as stated in this section) patent license to make, have made, - use, offer to sell, sell, import, and otherwise transfer the Work, - where such license applies only to those patent claims licensable - by such Contributor that are necessarily infringed by their - Contribution(s) alone or by combination of their Contribution(s) - with the Work to which such Contribution(s) was submitted. If You - institute patent litigation against any entity (including a - cross-claim or counterclaim in a lawsuit) alleging that the Work - or a Contribution incorporated within the Work constitutes direct - or contributory patent infringement, then any patent licenses - granted to You under this License for that Work shall terminate - as of the date such litigation is filed. - - 4. Redistribution. You may reproduce and distribute copies of the - Work or Derivative Works thereof in any medium, with or without - modifications, and in Source or Object form, provided that You - meet the following conditions: - - (a) You must give any other recipients of the Work or - Derivative Works a copy of this License; and - - (b) You must cause any modified files to carry prominent notices - stating that You changed the files; and - - (c) You must retain, in the Source form of any Derivative Works - that You distribute, all copyright, patent, trademark, and - attribution notices from the Source form of the Work, - excluding those notices that do not pertain to any part of - the Derivative Works; and - - (d) If the Work includes a "NOTICE" text file as part of its - distribution, then any Derivative Works that You distribute must - include a readable copy of the attribution notices contained - within such NOTICE file, excluding those notices that do not - pertain to any part of the Derivative Works, in at least one - of the following places: within a NOTICE text file distributed - as part of the Derivative Works; within the Source form or - documentation, if provided along with the Derivative Works; or, - within a display generated by the Derivative Works, if and - wherever such third-party notices normally appear. The contents - of the NOTICE file are for informational purposes only and - do not modify the License. You may add Your own attribution - notices within Derivative Works that You distribute, alongside - or as an addendum to the NOTICE text from the Work, provided - that such additional attribution notices cannot be construed - as modifying the License. - - You may add Your own copyright statement to Your modifications and - may provide additional or different license terms and conditions - for use, reproduction, or distribution of Your modifications, or - for any such Derivative Works as a whole, provided Your use, - reproduction, and distribution of the Work otherwise complies with - the conditions stated in this License. - - 5. Submission of Contributions. Unless You explicitly state otherwise, - any Contribution intentionally submitted for inclusion in the Work - by You to the Licensor shall be under the terms and conditions of - this License, without any additional terms or conditions. - Notwithstanding the above, nothing herein shall supersede or modify - the terms of any separate license agreement you may have executed - with Licensor regarding such Contributions. - - 6. Trademarks. This License does not grant permission to use the trade - names, trademarks, service marks, or product names of the Licensor, - except as required for reasonable and customary use in describing the - origin of the Work and reproducing the content of the NOTICE file. - - 7. Disclaimer of Warranty. Unless required by applicable law or - agreed to in writing, Licensor provides the Work (and each - Contributor provides its Contributions) on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or - implied, including, without limitation, any warranties or conditions - of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A - PARTICULAR PURPOSE. You are solely responsible for determining the - appropriateness of using or redistributing the Work and assume any - risks associated with Your exercise of permissions under this License. - - 8. Limitation of Liability. In no event and under no legal theory, - whether in tort (including negligence), contract, or otherwise, - unless required by applicable law (such as deliberate and grossly - negligent acts) or agreed to in writing, shall any Contributor be - liable to You for damages, including any direct, indirect, special, - incidental, or consequential damages of any character arising as a - result of this License or out of the use or inability to use the - Work (including but not limited to damages for loss of goodwill, - work stoppage, computer failure or malfunction, or any and all - other commercial damages or losses), even if such Contributor - has been advised of the possibility of such damages. - - 9. Accepting Warranty or Additional Liability. While redistributing - the Work or Derivative Works thereof, You may choose to offer, - and charge a fee for, acceptance of support, warranty, indemnity, - or other liability obligations and/or rights consistent with this - License. However, in accepting such obligations, You may act only - on Your own behalf and on Your sole responsibility, not on behalf - of any other Contributor, and only if You agree to indemnify, - defend, and hold each Contributor harmless for any liability - incurred by, or claims asserted against, such Contributor by reason - of your accepting any such warranty or additional liability. - - END OF TERMS AND CONDITIONS - - APPENDIX: How to apply the Apache License to your work. - - To apply the Apache License to your work, attach the following - boilerplate notice, with the fields enclosed by brackets "[]" - replaced with your own identifying information. (Don't include - the brackets!) The text should be enclosed in the appropriate - comment syntax for the file format. We also recommend that a - file or class name and description of purpose be included on the - same "printed page" as the copyright notice for easier - identification within third-party archives. - - Copyright [yyyy] [name of copyright owner] - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. diff --git a/ThirdParty/ITK/N4 b/ThirdParty/ITK/N4 deleted file mode 100755 index 2839d98..0000000 Binary files a/ThirdParty/ITK/N4 and /dev/null differ diff --git a/pipelines/neonatal-pipeline-v1.sh b/pipelines/neonatal-pipeline-v1.sh deleted file mode 100755 index 0ef1d7e..0000000 --- a/pipelines/neonatal-pipeline-v1.sh +++ /dev/null @@ -1,145 +0,0 @@ -#! /bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - -usage() -{ - base=$(basename "$0") - echo "usage: $base subjectID [options] -This script runs the neonatal segmentation pipeline of Draw-EM. -Required Options: - -a / -age Subject age in weeks. This is used to select the appropriate template for the initial registration. - If the age is <28w or >44w, it will be set to 28w or 44w respectively. - - Note: The age of the subject can be alternatively provided by a file named 'ages.csv' inside the data directory. - The format of the file must be a space-delimited file with lines: 'subjectID age' - -Additional options: - -c / -cleanup <0/1> Whether cleanup of temporary files is required (default: 1) - -d / -data-dir The directory used to run the script. - It must contain a T2 folder and the T2/subjectID.nii.gz image of the subject. - -p / -save-posteriors <0/1> Whether the structures' posteriors are required (default: 0) - -t / -threads Number of threads (CPU cores) allowed for the registration to run in parallel (default: 1) - -v / -verbose <0/1> Whether the script progress is reported (default: 1) - -h / -help / --help Print usage. -" - exit; -} - - - -if [ -n "$DRAWEMDIR" ]; then - [ -d "$DRAWEMDIR" ] || { echo "DRAWEMDIR environment variable invalid!" 1>&2; exit 1; } -else - export DRAWEMDIR="$(cd "$(dirname "$BASH_SOURCE")"/.. && pwd)" -fi - -[ $# -ge 1 ] || { usage; } -subj=$1 - -case "$subj" in - -h|-help|--help) usage; ;; -esac - -age="" -cleanup=1 # whether to delete temporary files once done -datadir=`pwd` -posteriors=0 # whether to output posterior probability maps -threads=1 -verbose=1 - -atlasname=non-rigid-v2 -version="v1" - -while [ $# -gt 0 ]; do - case "$2" in - -a|-age) shift; age=$2; ;; - -c|-cleanup) shift; cleanup=$2; ;; - -d|-data-dir) shift; datadir=$2; ;; - -p|-save-posteriors) shift; posteriors=$2; ;; - -t|-threads) shift; threads=$2; ;; - -v|-verbose) shift; verbose=$2; ;; - -h|-help|--help) usage; ;; - -*) echo "$0: Unrecognized option $1" >&2; usage; ;; - *) break ;; - esac - shift -done - - -cd "$datadir" || { echo "$0: Directory $datadir does not exist" >&2;exit 1; } - -[ -n "$age" -o ! -f ages.csv ] || { age=`cat ages.csv | tr '\t' ' ' | grep "^$subj " | cut -d' ' -f2`; } -[ -n "$age" ] || { echo "Subject age must be provided as argument or from the $datadir/ages.csv file" >&2; exit 1; } -age=`printf "%.*f\n" 0 $age` #round -[ $age -gt 28 ] || { age=28; } -[ $age -lt 44 ] || { age=44; } - - -[ $verbose -le 0 ] || { echo "DrawEM multi atlas $version -Subject: $subj -Age: $age -Directory: $datadir -Posteriors: $posteriors -Cleanup: $cleanup -Threads: $threads - -$BASH_SOURCE $@ -----------------------------"; } - -mkdir -p logs || exit 1 - -run() -{ - [ $verbose -le 0 ] || echo -n "$@..." - "$DRAWEMDIR/scripts/$version/$@" 1>>logs/$subj 2>>logs/$subj-err - if [ $? -eq 0 ]; then - [ $verbose -le 0 ] || echo " done" - else - [ $verbose -le 0 ] || echo " failed: see log file logs/$subj-err for details" - exit 1 - fi -} - -rm -f logs/$subj logs/$subj-err -run preprocess.sh $subj $age -# phase 1 tissue segmentation -run tissue-priors.sh $subj $age $atlasname $threads -# registration using gm posterior + image -run register-multi-atlas-using-gm-posteriors.sh $subj $age $threads -# structural segmentation -run labels-multi-atlas.sh $subj -run segmentation.sh $subj -# post-processing -run separate-hemispheres.sh $subj -run correct-segmentation.sh $subj -run postprocess.sh $subj - -# if probability maps are required -[ "$posteriors" == "0" -o "$posteriors" == "no" -o "$posteriors" == "false" ] || run postprocess-pmaps.sh $subj - -# cleanup -if [ "$cleanup" == "1" -o "$cleanup" == "yes" -o "$cleanup" == "true" ] && [ -f "segmentations/${subj}_labels.nii.gz" ];then - run clear-data.sh $subj - rm -f logs/$subj logs/$subj-err - rmdir logs 2> /dev/null # may fail if other log files exist -fi - -exit 0 diff --git a/scripts/v1.1/clear-data.sh b/scripts/clear-data.sh similarity index 100% rename from scripts/v1.1/clear-data.sh rename to scripts/clear-data.sh diff --git a/scripts/v1.1/clear-small-components.sh b/scripts/clear-small-components.sh similarity index 98% rename from scripts/v1.1/clear-small-components.sh rename to scripts/clear-small-components.sh index cf089ba..581267e 100755 --- a/scripts/v1.1/clear-small-components.sh +++ b/scripts/clear-small-components.sh @@ -24,11 +24,6 @@ outf=$2 keepratio=0.05 if [ $# -gt 2 ];then keepratio=$3;fi -run(){ - echo "$@" - "$@" || exit 1 -} - if [ ! -f $outf ];then sdir=segmentations-data diff --git a/scripts/v1.1/clear-small-ventricle-components-labels.sh b/scripts/clear-small-ventricle-components-labels.sh similarity index 98% rename from scripts/v1.1/clear-small-ventricle-components-labels.sh rename to scripts/clear-small-ventricle-components-labels.sh index d1f4c84..c7f98e1 100755 --- a/scripts/v1.1/clear-small-ventricle-components-labels.sh +++ b/scripts/clear-small-ventricle-components-labels.sh @@ -26,11 +26,6 @@ scriptdir=$(dirname "$BASH_SOURCE") sdir=segmentations-data mkdir -p $sdir/corrections || exit 1 -run(){ - echo "$@" - "$@" || exit 1 -} - # cleaning up small ventricle components run mirtk padding segmentations/$subj-em.nii.gz segmentations/$subj-em.nii.gz $sdir/corrections/$subj-hwm-init.nii.gz 26 0 -invert run mirtk padding segmentations/$subj-initial.nii.gz segmentations/$subj-initial.nii.gz $sdir/corrections/$subj-ven-init.nii.gz 2 49 50 0 -invert 2 49 50 1 diff --git a/scripts/v1/correct-hemisphere-holes.sh b/scripts/correct-hemisphere-holes.sh similarity index 98% rename from scripts/v1/correct-hemisphere-holes.sh rename to scripts/correct-hemisphere-holes.sh index 9b8c5e0..71c856c 100755 --- a/scripts/v1/correct-hemisphere-holes.sh +++ b/scripts/correct-hemisphere-holes.sh @@ -27,11 +27,6 @@ mkdir -p $sdir/corrections || exit 1 subcorts=`cat $DRAWEMDIR/parameters/subcortical-all.csv` -run(){ - echo "$@" - "$@" || exit 1 -} - if [ -f segmentations/${subj}_L_white.nii.gz -a -f segmentations/${subj}_R_white.nii.gz ];then # if small gm components exist inside the white surface.. diff --git a/scripts/v1.1/correct-segmentation.sh b/scripts/correct-segmentation.sh similarity index 100% rename from scripts/v1.1/correct-segmentation.sh rename to scripts/correct-segmentation.sh diff --git a/scripts/v1.1/labels-multi-atlas.sh b/scripts/labels-multi-atlas.sh similarity index 99% rename from scripts/v1.1/labels-multi-atlas.sh rename to scripts/labels-multi-atlas.sh index ca579f5..154edb2 100755 --- a/scripts/v1.1/labels-multi-atlas.sh +++ b/scripts/labels-multi-atlas.sh @@ -51,7 +51,7 @@ atlases="" num=0 #for each atlas -for atnum in {01..20};do +for atnum in 0{1..9} {10..20};do atlas="ALBERT_"$atnum if [ ! -f dofs/$subj-$atlas-n.dof.gz ];then continue;fi diff --git a/scripts/v1.1/postprocess-cortical.sh b/scripts/postprocess-cortical.sh similarity index 100% rename from scripts/v1.1/postprocess-cortical.sh rename to scripts/postprocess-cortical.sh diff --git a/scripts/v1/postprocess-pmaps.sh b/scripts/postprocess-pmaps.sh similarity index 95% rename from scripts/v1/postprocess-pmaps.sh rename to scripts/postprocess-pmaps.sh index a3886f2..afbd628 100755 --- a/scripts/v1/postprocess-pmaps.sh +++ b/scripts/postprocess-pmaps.sh @@ -23,12 +23,6 @@ [ $# -eq 1 ] || { echo "usage: $(basename "$0") " 1>&2; exit 1; } subj=$1 - -run(){ - echo "$@" - "$@" || exit 1 -} - rdir=posteriors sdir=segmentations-data @@ -49,6 +43,11 @@ for ((n=0;n<${#all[*]};n++));do r=${all[$n]}; mkdir -p $rdir/seg$r || exit 1;don cp $sdir/posteriors/csf/$subj.nii.gz $rdir/seg83/$subj.nii.gz || exit 1 cp $sdir/posteriors/outlier/$subj.nii.gz $rdir/seg84/$subj.nii.gz || exit 1 +for tiss in csf gm wm;do + mkdir -p $rdir/$tiss + cp $sdir/posteriors/gm/$subj.nii.gz $rdir/$tiss/$subj.nii.gz || exit 1 +done + # cortical wm, gm addem="" for ((n=0;n<$numcorts;n++));do r=${corts[$n]}; addem=$addem"-add $sdir/labels/seg$r-extended/$subj.nii.gz "; done diff --git a/scripts/v1/postprocess.sh b/scripts/postprocess.sh similarity index 98% rename from scripts/v1/postprocess.sh rename to scripts/postprocess.sh index b2ef58c..763e920 100755 --- a/scripts/v1/postprocess.sh +++ b/scripts/postprocess.sh @@ -32,11 +32,6 @@ scriptdir=$(dirname "$BASH_SOURCE") f=segmentations/$subj-initial.nii.gz -run(){ - echo "$@" - "$@" || exit 1 -} - if [ ! -f segmentations/"$subj"_all_labels$suffix.nii.gz ];then # creating the all labels file (initial segmentation + cortical division) $scriptdir/postprocess-cortical.sh $subj diff --git a/scripts/v1.1/preprocess.sh b/scripts/preprocess.sh similarity index 88% rename from scripts/v1.1/preprocess.sh rename to scripts/preprocess.sh index 8093709..ed3f467 100755 --- a/scripts/v1.1/preprocess.sh +++ b/scripts/preprocess.sh @@ -32,12 +32,6 @@ else export PATH="$FSLDIR/bin:$PATH" fi -run(){ - echo "$@" - "$@" || exit 1 -} - - sdir=segmentations-data mkdir -p segmentations N4 dofs bias || exit 1 @@ -53,7 +47,7 @@ if [ ! -f N4/$subj.nii.gz ];then fi #bias correct - run $DRAWEMDIR/ThirdParty/ITK/N4 3 -i N4/${subj}_rescaled.nii.gz -x segmentations/${subj}_brain_mask.nii.gz -o "[N4/${subj}_corr.nii.gz,bias/$subj.nii.gz]" -c "[50x50x50,0.001]" -s 2 -b "[100,3]" -t "[0.15,0.01,200]" + run mirtk N4 3 -i N4/${subj}_rescaled.nii.gz -x segmentations/${subj}_brain_mask.nii.gz -o "[N4/${subj}_corr.nii.gz,bias/$subj.nii.gz]" -c "[50x50x50,0.001]" -s 2 -b "[100,3]" -t "[0.15,0.01,200]" run mirtk calculate N4/${subj}_corr.nii.gz -mul segmentations/${subj}_brain_mask.nii.gz -out N4/${subj}_corr.nii.gz #rescale image diff --git a/scripts/v1.1/register-multi-atlas-using-gm-posteriors.sh b/scripts/register-multi-atlas-using-gm-posteriors.sh similarity index 93% rename from scripts/v1.1/register-multi-atlas-using-gm-posteriors.sh rename to scripts/register-multi-atlas-using-gm-posteriors.sh index 19f0218..1272f20 100755 --- a/scripts/v1.1/register-multi-atlas-using-gm-posteriors.sh +++ b/scripts/register-multi-atlas-using-gm-posteriors.sh @@ -26,21 +26,15 @@ age=$2 njobs=1 if [ $# -gt 2 ];then njobs=$3;fi -run(){ - echo "$@" - "$@" || exit 1 -} - -scriptdir=$(dirname "$BASH_SOURCE") sdir=segmentations-data mkdir -p dofs -for i in {01..20};do +for i in 0{1..9} {10..20};do atlas=ALBERT_$i dof=dofs/$subj-$atlas-n.dof.gz if [ ! -f $dof ];then - run mirtk register N4/$subj.nii.gz $DRAWEMDIR/atlases/ALBERTs/T2/$atlas.nii.gz $sdir/gm-posteriors/$subj.nii.gz $DRAWEMDIR/atlases/ALBERTs/gm-posteriors-v3/$atlas.nii.gz -parin $DRAWEMDIR/parameters/ireg-multichannel-structural.cfg -dofout $dof -threads $njobs + run mirtk register N4/$subj.nii.gz $DRAWEMDIR/atlases/ALBERTs/T2/$atlas.nii.gz $sdir/gm-posteriors/$subj.nii.gz $DRAWEMDIR/atlases/ALBERTs/gm-posteriors-v3/$atlas.nii.gz -parin $DRAWEMDIR/parameters/ireg-multichannel-structural.cfg -dofout $dof -threads $njobs -v 0 fi done diff --git a/scripts/v1.1/segmentation.sh b/scripts/segmentation.sh similarity index 98% rename from scripts/v1.1/segmentation.sh rename to scripts/segmentation.sh index fb04f87..e959965 100755 --- a/scripts/v1.1/segmentation.sh +++ b/scripts/segmentation.sh @@ -23,11 +23,6 @@ [ $# -eq 1 ] || { echo "usage: $(basename "$0") "; exit 1; } subj=$1 -run() -{ - echo "$@" - "$@" || exit 1 -} if [ ! -f segmentations/$subj-initial.nii.gz ];then sdir=segmentations-data diff --git a/scripts/v1.1/separate-hemispheres.sh b/scripts/separate-hemispheres.sh similarity index 99% rename from scripts/v1.1/separate-hemispheres.sh rename to scripts/separate-hemispheres.sh index 7c2ebdd..d1f058e 100755 --- a/scripts/v1.1/separate-hemispheres.sh +++ b/scripts/separate-hemispheres.sh @@ -22,11 +22,6 @@ [ $# -eq 1 ] || { echo "usage: $(basename "$0") "; exit 1; } subj=$1 -run(){ - echo "$@" - "$@" || exit 1 -} - scriptdir=$(dirname "$BASH_SOURCE") diff --git a/scripts/v1.1/tissue-priors.sh b/scripts/tissue-priors.sh similarity index 95% rename from scripts/v1.1/tissue-priors.sh rename to scripts/tissue-priors.sh index 2c1b7e5..fee8daa 100755 --- a/scripts/v1.1/tissue-priors.sh +++ b/scripts/tissue-priors.sh @@ -27,22 +27,16 @@ atlasname=$3 njobs=1 if [ $# -gt 3 ];then njobs=$4;fi -run(){ - echo "$@" - "$@" || exit 1 -} - - #registration dof=dofs/$subj-template-$age-n.dof.gz -mkdir -p $(dirname "${dof}") +mkdir -p dofs if [ ! -f $dof ];then - run mirtk register N4/$subj.nii.gz $DRAWEMDIR/atlases/$atlasname/T2/template-$age.nii.gz -dofout $dof -parin $DRAWEMDIR/parameters/ireg.cfg -threads $njobs + run mirtk register N4/$subj.nii.gz $DRAWEMDIR/atlases/$atlasname/T2/template-$age.nii.gz -dofout $dof -parin $DRAWEMDIR/parameters/ireg.cfg -threads $njobs -v 0 fi needtseg=0 -for i in {01..20};do +for i in 0{1..9} {10..20};do atlas=ALBERT_$i if [ ! -f dofs/$subj-$atlas-n.dof.gz ];then needtseg=1 diff --git a/scripts/v1.1/correct-hemisphere-holes.sh b/scripts/v1.1/correct-hemisphere-holes.sh deleted file mode 100755 index 9b8c5e0..0000000 --- a/scripts/v1.1/correct-hemisphere-holes.sh +++ /dev/null @@ -1,63 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $BASH_SOURCE " 1>&2; exit 1; } -subj=$1 - -sdir=segmentations-data -mkdir -p $sdir/corrections || exit 1 - -subcorts=`cat $DRAWEMDIR/parameters/subcortical-all.csv` - -run(){ - echo "$@" - "$@" || exit 1 -} - -if [ -f segmentations/${subj}_L_white.nii.gz -a -f segmentations/${subj}_R_white.nii.gz ];then - -# if small gm components exist inside the white surface.. -run mirtk calculate segmentations/$subj-initial.nii.gz -mul 0 -out $sdir/corrections/$subj-gmtochange.nii.gz -for h in L R;do -run mirtk calculate segmentations/${subj}_${h}_white.nii.gz -mul segmentations/$subj-initial.nii.gz -out $sdir/corrections/$subj-gmmask.nii.gz -run mirtk padding $sdir/corrections/$subj-gmtochange.nii.gz $sdir/corrections/$subj-gmmask.nii.gz $sdir/corrections/$subj-gmtochange.nii.gz 1000 1 -done -rm $sdir/corrections/$subj-gmmask.nii.gz - -volcorr=`mirtk measure-volume $sdir/corrections/$subj-gmtochange.nii.gz` -if [ "$volcorr" != "" ];then - -# ..redistribute small components' probability into wm/dgm -num=1; -inprobs="$sdir/posteriors/wm/$subj.nii.gz 2000" -outprobs="$sdir/posteriors/wm/$subj.nii.gz" -for r in ${subcorts}; do -inprobs="$inprobs $sdir/posteriors/seg$r/$subj.nii.gz $r" -outprobs="$outprobs $sdir/posteriors/seg$r/$subj.nii.gz" -let num=num+1 -done - -run mirtk change-label segmentations/$subj-initial.nii.gz $sdir/corrections/$subj-gmtochange.nii.gz $num $inprobs segmentations/$subj-initial.nii.gz $sdir/posteriors/gm/$subj.nii.gz $outprobs -run mirtk padding $sdir/posteriors/gm/$subj.nii.gz $sdir/corrections/$subj-gmtochange.nii.gz $sdir/posteriors/gm/$subj.nii.gz 1 0 - -fi - -fi diff --git a/scripts/v1.1/postprocess-pmaps.sh b/scripts/v1.1/postprocess-pmaps.sh deleted file mode 100755 index a3886f2..0000000 --- a/scripts/v1.1/postprocess-pmaps.sh +++ /dev/null @@ -1,70 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") " 1>&2; exit 1; } -subj=$1 - - -run(){ - echo "$@" - "$@" || exit 1 -} - -rdir=posteriors -sdir=segmentations-data - - -subcorts=(`cat $DRAWEMDIR/parameters/subcortical-all.csv`) -corts=(`cat $DRAWEMDIR/parameters/cortical.csv`) -wmcorts=(`cat $DRAWEMDIR/parameters/cortical-wm.csv`) -all=(`cat $DRAWEMDIR/parameters/all-labels.csv `) -numcorts=${#corts[*]} -numsubcorts=${#subcorts[*]} - - -mkdir -p $sdir/posteriors/temp || exit 1 -for ((n=0;n<${#all[*]};n++));do r=${all[$n]}; mkdir -p $rdir/seg$r || exit 1;done - - -# out, csf -cp $sdir/posteriors/csf/$subj.nii.gz $rdir/seg83/$subj.nii.gz || exit 1 -cp $sdir/posteriors/outlier/$subj.nii.gz $rdir/seg84/$subj.nii.gz || exit 1 - -# cortical wm, gm -addem="" -for ((n=0;n<$numcorts;n++));do r=${corts[$n]}; addem=$addem"-add $sdir/labels/seg$r-extended/$subj.nii.gz "; done -addem=`echo $addem|sed -e 's:^-add::g'` -run mirtk calculate $addem -out $sdir/posteriors/temp/$subj-gmwm.nii.gz - -for ((n=0;n<$numcorts;n++));do -r=${corts[$n]}; -val=${wmcorts[$n]}; -run mirtk calculate $sdir/labels/seg$r-extended/$subj.nii.gz -div-with-zero $sdir/posteriors/temp/$subj-gmwm.nii.gz -mul $sdir/posteriors/gm/$subj.nii.gz -out $rdir/seg$r/$subj.nii.gz -run mirtk calculate $sdir/labels/seg$r-extended/$subj.nii.gz -div-with-zero $sdir/posteriors/temp/$subj-gmwm.nii.gz -mul $sdir/posteriors/wm/$subj.nii.gz -out $rdir/seg$val/$subj.nii.gz -done -rm $sdir/posteriors/temp/$subj-gmwm.nii.gz - -# subcortical -for ((n=0;n<$numsubcorts;n++));do -r=${subcorts[$n]}; -cp $sdir/posteriors/seg$r/$subj.nii.gz $rdir/seg$r/$subj.nii.gz || exit 1 -done diff --git a/scripts/v1.1/postprocess.sh b/scripts/v1.1/postprocess.sh deleted file mode 100755 index b2ef58c..0000000 --- a/scripts/v1.1/postprocess.sh +++ /dev/null @@ -1,68 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -ge 1 -a $# -le 2 ] || { echo "usage: $(basename "$0") []" 1>&2; exit 1; } - -subj=$1 -suffix="" -[ $# -eq 1 ] || suffix=$2 - -sdir=segmentations-data -scriptdir=$(dirname "$BASH_SOURCE") - -f=segmentations/$subj-initial.nii.gz - - -run(){ - echo "$@" - "$@" || exit 1 -} - -if [ ! -f segmentations/"$subj"_all_labels$suffix.nii.gz ];then -# creating the all labels file (initial segmentation + cortical division) -$scriptdir/postprocess-cortical.sh $subj -run mirtk padding $sdir/cortical/$subj.nii.gz $f segmentations/$subj-cortical-wm.nii.gz 2000 0 -invert -run mirtk padding $sdir/cortical/$subj.nii.gz $f segmentations/$subj-cortical-gm.nii.gz 1000 0 -invert - -corts=(`cat $DRAWEMDIR/parameters/cortical.csv`) -wmcorts=(`cat $DRAWEMDIR/parameters/cortical-wm.csv`) -gmtowmnumbers="" -for ((r=0;r<${#corts[*]};r++));do -gmtowmnumbers="$gmtowmnumbers 1 ${corts[$r]} ${wmcorts[$r]}" -done -run mirtk padding segmentations/$subj-cortical-wm.nii.gz segmentations/$subj-cortical-wm.nii.gz segmentations/$subj-cortical-wm.nii.gz $gmtowmnumbers -run mirtk padding $f $f segmentations/"$subj"_all_labels_ini$suffix.nii.gz 2 1000 2000 0 -run mirtk calculate segmentations/"$subj"_all_labels_ini$suffix.nii.gz -add segmentations/$subj-cortical-gm.nii.gz -add segmentations/$subj-cortical-wm.nii.gz -out segmentations/"$subj"_all_labels$suffix.nii.gz -# cleanup -rm segmentations/$subj-cortical-gm.nii.gz segmentations/$subj-cortical-wm.nii.gz segmentations/"$subj"_all_labels_ini$suffix.nii.gz -fi - -if [ ! -f segmentations/"$subj"_labels$suffix.nii.gz ];then -# creating the labels file -run mirtk padding segmentations/"$subj"_all_labels$suffix.nii.gz segmentations/"$subj"_all_labels$suffix.nii.gz segmentations/"$subj"_labels$suffix.nii.gz $DRAWEMDIR/parameters/all-labels-to-labels.csv -fi - -if [ ! -f segmentations/"$subj"_tissue_labels$suffix.nii.gz ];then -# creating the tissue labels labels -run mirtk padding segmentations/"$subj"_all_labels$suffix.nii.gz segmentations/"$subj"_all_labels$suffix.nii.gz segmentations/"$subj"_tissue_labels$suffix.nii.gz $DRAWEMDIR/parameters/all-labels-to-tissue-labels.csv -fi - diff --git a/scripts/v1/clear-data.sh b/scripts/v1/clear-data.sh deleted file mode 100755 index d08af34..0000000 --- a/scripts/v1/clear-data.sh +++ /dev/null @@ -1,74 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") " 1>&2; exit 1; } -subj=$1 - -sdir=segmentations-data -[ -n "$sdir" ] || { echo "$BASH_SOURCE: sdir must not be empty!" 1>&2; exit 1; } - -rm -f $sdir/MADs/$subj.nii.gz $sdir/MADs/$subj-grad.nii.gz $sdir/MADs/$subj-subspace.nii.gz -rm -f $sdir/atlas-weights/$subj-ALBERT_*.nii.gz -rm -f $sdir/brain/$subj.nii.gz $sdir/brain/${subj}_brain.nii.gz $sdir/brain/${subj}_brain_mask.nii.gz -rm -f $sdir/corrections/$subj-gmtochange.nii.gz $sdir/corrections/$subj-ventohwm.nii.gz -rm -f $sdir/cortical/$subj.nii.gz -rm -f $sdir/gm-posteriors/$subj.nii.gz -rm -f $sdir/labels/*/$subj.nii.gz -rm -f $sdir/posteriors/*/$subj.nii.gz -rm -f $sdir/template/*/$subj.nii.gz -rm -f $sdir/tissue-initial-segmentations/$subj.nii.gz -rm -f $sdir/tissue-posteriors/*/$subj.nii.gz -rm -f $sdir/transformations/T2-$subj-ALBERT_*.nii.gz $sdir/transformations/$subj-ALBERT_*.nii.gz -rm -f segmentations/$subj-em.nii.gz segmentations/$subj-initial.nii.gz -rm -f logs/$subj logs/$subj-err logs/$subj-em logs/$subj-em-err logs/$subj-tissue-em logs/$subj-tissue-em-err - -#may fail if other files present, which is ok -rmdir $sdir/atlas-weights 2> /dev/null -rmdir $sdir/brain 2> /dev/null -rmdir $sdir/corrections 2> /dev/null -rmdir $sdir/cortical 2> /dev/null -rmdir $sdir/gm-posteriors 2> /dev/null -rmdir $sdir/labels/* 2> /dev/null -rmdir $sdir/labels 2> /dev/null -rmdir $sdir/MADs 2> /dev/null -rmdir $sdir/posteriors/* 2> /dev/null -rmdir $sdir/posteriors 2> /dev/null -rmdir $sdir/template/* 2> /dev/null -rmdir $sdir/template 2> /dev/null -rmdir $sdir/tissue-initial-segmentations 2> /dev/null -rmdir $sdir/tissue-posteriors/* 2> /dev/null -rmdir $sdir/tissue-posteriors 2> /dev/null -rmdir $sdir/transformations 2> /dev/null -rmdir $sdir 2> /dev/null -rmdir logs 2> /dev/null - -exit 0 - - - - - - - - - - - diff --git a/scripts/v1/clear-small-components.sh b/scripts/v1/clear-small-components.sh deleted file mode 100755 index cf089ba..0000000 --- a/scripts/v1/clear-small-components.sh +++ /dev/null @@ -1,58 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - -[ $# -ge 2 ] || { echo "usage: $BASH_SOURCE []" 1>&2; exit 1; } -f=$1 -outf=$2 -keepratio=0.05 -if [ $# -gt 2 ];then keepratio=$3;fi - -run(){ - echo "$@" - "$@" || exit 1 -} - -if [ ! -f $outf ];then - -sdir=segmentations-data -mkdir -p $sdir/corrections || exit 1 -base=`basename $f| sed -e 's:.nii.gz::g'` - -# connected components -run mirtk extract-connected-components $f $sdir/corrections/$base-lccs.nii.gz -all -output-component-labels -connectivity 6 -mirtk measure-volume $sdir/corrections/$base-lccs.nii.gz > $sdir/corrections/$base-lccs-comps -comp1=`cat $sdir/corrections/$base-lccs-comps |grep "^1 "|cut -d' ' -f2` -comps=`cat $sdir/corrections/$base-lccs-comps |wc -l` -retainvol=`echo "$comp1*$keepratio" | /usr/bin/bc` -retainvol=`printf "%.*f\n" 0 $retainvol` #round -retain="1" - -# keep only connected components with volume > 5% volume of first component -for ((r=2;r<$comps;r++));do - comp=`cat $sdir/corrections/$base-lccs-comps |grep "^$r "|cut -d' ' -f2` - if [ `echo "$comp<$retainvol" | /usr/bin/bc` -eq 1 ];then break;fi - retain="$retain $r" -done - -numretain=`echo $retain|wc -w` -run mirtk padding $f $sdir/corrections/$base-lccs.nii.gz $outf $numretain $retain 0 -invert - -rm $sdir/corrections/$base-lccs.nii.gz $sdir/corrections/$base-lccs-comps -fi diff --git a/scripts/v1/clear-small-ventricle-components-labels.sh b/scripts/v1/clear-small-ventricle-components-labels.sh deleted file mode 100755 index ad60cba..0000000 --- a/scripts/v1/clear-small-ventricle-components-labels.sh +++ /dev/null @@ -1,63 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -ge 1 ] || { echo "usage: $BASH_SOURCE " 1>&2; exit 1; } -subj=$1 - -scriptdir=$(dirname "$BASH_SOURCE") -sdir=segmentations-data -mkdir -p $sdir/corrections || exit 1 - -run(){ - echo "$@" - "$@" || exit 1 -} - -# cleaning up small ventricle components -run mirtk padding segmentations/$subj-em.nii.gz segmentations/$subj-em.nii.gz $sdir/corrections/$subj-hwm-init.nii.gz 26 0 -invert -run mirtk padding segmentations/$subj-initial.nii.gz segmentations/$subj-initial.nii.gz $sdir/corrections/$subj-ven-init.nii.gz 2 49 50 0 -invert 2 49 50 1 -$scriptdir/clear-small-components.sh $sdir/corrections/$subj-ven-init.nii.gz $sdir/corrections/$subj-ven.nii.gz -run mirtk calculate $sdir/corrections/$subj-ven-init.nii.gz -sub $sdir/corrections/$subj-ven.nii.gz -out $sdir/corrections/$subj-ven-diff.nii.gz - -# if small ventricle components are surrounded by hwm.. -run mirtk fill-holes $sdir/corrections/$subj-hwm-init.nii.gz $sdir/corrections/$subj-hwm-fillh.nii.gz -run mirtk calculate $sdir/corrections/$subj-ven-diff.nii.gz -mul $sdir/corrections/$subj-hwm-fillh.nii.gz -out $sdir/corrections/$subj-ventohwm.nii.gz - -volcorr=`mirtk measure-volume $sdir/corrections/$subj-ventohwm.nii.gz` -if [ "$volcorr" != "" ];then -# ..redistribute small ventricle components' probability into hwm/wm -run mirtk calculate $sdir/posteriors/seg49/$subj.nii.gz -add $sdir/posteriors/seg50/$subj.nii.gz -mul $sdir/corrections/$subj-ventohwm.nii.gz -out $sdir/corrections/$subj-ventohwm-prob.nii.gz -run mirtk calculate $sdir/corrections/$subj-ventohwm-prob.nii.gz -add $sdir/posteriors/hwm/$subj.nii.gz -out $sdir/posteriors/hwm/$subj.nii.gz -run mirtk calculate $sdir/corrections/$subj-ventohwm-prob.nii.gz -add $sdir/posteriors/wm/$subj.nii.gz -out $sdir/posteriors/wm/$subj.nii.gz -run mirtk padding $sdir/posteriors/seg49/$subj.nii.gz $sdir/corrections/$subj-ventohwm.nii.gz $sdir/posteriors/seg49/$subj.nii.gz 1 0 -run mirtk padding $sdir/posteriors/seg50/$subj.nii.gz $sdir/corrections/$subj-ventohwm.nii.gz $sdir/posteriors/seg50/$subj.nii.gz 1 0 - -# ..fix the segmentation too -run mirtk padding segmentations/$subj-initial.nii.gz $sdir/corrections/$subj-ventohwm.nii.gz segmentations/$subj-initial.nii.gz 1 2000 - -# clean up -rm $sdir/corrections/$subj-ventohwm-prob.nii.gz -fi - -# clean up -rm $sdir/corrections/$subj-ven.nii.gz $sdir/corrections/$subj-ven-diff.nii.gz $sdir/corrections/$subj-hwm-init.nii.gz $sdir/corrections/$subj-hwm-fillh.nii.gz $sdir/corrections/$subj-ven-init.nii.gz - - diff --git a/scripts/v1/correct-segmentation.sh b/scripts/v1/correct-segmentation.sh deleted file mode 100755 index 84c5fa4..0000000 --- a/scripts/v1/correct-segmentation.sh +++ /dev/null @@ -1,33 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") "; exit 1; } -subj=$1 - -scriptdir=$(dirname "$BASH_SOURCE") - -# correct for small ventricle components -$scriptdir/clear-small-ventricle-components-labels.sh $subj -# correct holes in the segmentation of the hemispheres -$scriptdir/correct-hemisphere-holes.sh $subj - - - diff --git a/scripts/v1/labels-multi-atlas.sh b/scripts/v1/labels-multi-atlas.sh deleted file mode 100755 index d26d086..0000000 --- a/scripts/v1/labels-multi-atlas.sh +++ /dev/null @@ -1,153 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") " 1>&2; exit 1; } -subj=$1 - -sdir=segmentations-data - -corts=`cat $DRAWEMDIR/parameters/cortical.csv` -wmcorts=`cat $DRAWEMDIR/parameters/cortical-wm.csv` -subcorts=`cat $DRAWEMDIR/parameters/subcortical-all.csv` -all=`cat $DRAWEMDIR/parameters/all-labels.csv ` - - -run(){ - echo "$@" - "$@" || exit 1 -} - -if [ ! -f $sdir/MADs/$subj-subspace.nii.gz ];then - -mkdir -p $sdir/MADs $sdir/cortical $sdir/transformations $sdir/atlas-weights || exit 1 -for r in ${all};do mkdir -p $sdir/labels/seg$r || exit 1; done -sigma=10000 -run mirtk convert-image N4/$subj.nii.gz $sdir/atlas-weights/$subj-normalized.nii.gz -rescale 0 200 -double - - -atlases="" -num=0 - -#for each atlas -for atnum in {01..20};do -atlas="ALBERT_"$atnum -if [ ! -f dofs/$subj-$atlas-n.dof.gz ];then continue;fi - -#transform atlas labels -if [ ! -f $sdir/transformations/$subj-$atlas.nii.gz ];then -ms=$sdir/template/$str/$subj.nii.gz -run mirtk transform-image $DRAWEMDIR/atlases/ALBERTs/segmentations-v2.1/$atlas.nii.gz $sdir/transformations/$subj-$atlas.nii.gz -target N4/$subj.nii.gz -dofin dofs/$subj-$atlas-n.dof.gz -interp NN -fi - -#transform atlas labels -if [ ! -f $sdir/transformations/T2-$subj-$atlas.nii.gz ];then -ms=$sdir/template/$str/$subj.nii.gz -run mirtk transform-image $DRAWEMDIR/atlases/ALBERTs/T2/$atlas.nii.gz $sdir/transformations/T2-$subj-$atlas.nii.gz -target N4/$subj.nii.gz -dofin dofs/$subj-$atlas-n.dof.gz -interp BSpline -fi - -#weight atlas labels locally -if [ ! -f $sdir/atlas-weights/$subj-$atlas.nii.gz ];then -run mirtk normalize N4/$subj.nii.gz $sdir/transformations/T2-$subj-$atlas.nii.gz $sdir/atlas-weights/$subj-$atlas-normalized.nii.gz -piecewise -run mirtk convert-image $sdir/atlas-weights/$subj-$atlas-normalized.nii.gz $sdir/atlas-weights/$subj-$atlas-normalized.nii.gz -rescale 0 200 -double -# This calculates weights based on a gaussian distance. -run mirtk calculate $sdir/atlas-weights/$subj-normalized.nii.gz -sub $sdir/atlas-weights/$subj-$atlas-normalized.nii.gz -sq -out $sdir/atlas-weights/$subj-$atlas.nii.gz -run mirtk calculate-filtering $sdir/atlas-weights/$subj-$atlas.nii.gz -kernel 3 -mean $sdir/atlas-weights/$subj-$atlas.nii.gz -run mirtk calculate $sdir/atlas-weights/$subj-$atlas.nii.gz -mul -27 -div-with-zero $sigma -out $sdir/atlas-weights/$subj-$atlas.nii.gz -run mirtk calculate $sdir/atlas-weights/$subj-$atlas.nii.gz -exp -out $sdir/atlas-weights/$subj-$atlas.nii.gz -# smooth the weights -run mirtk calculate-filtering $sdir/atlas-weights/$subj-$atlas.nii.gz -kernel 3 -median $sdir/atlas-weights/$subj-$atlas.nii.gz -rm -f $sdir/atlas-weights/$subj-$atlas-normalized.nii.gz -fi - -atlases="$atlases $atlas" -num=$(($num+1)) -done - -rm -f $sdir/atlas-weights/$subj-normalized.nii.gz - - -#split labels -transformed="" -transformedw="" -for atlas in ${atlases};do -transformed="$transformed $sdir/transformations/$subj-$atlas.nii.gz"; -transformedw="$transformedw $sdir/atlas-weights/$subj-$atlas.nii.gz"; -done - -splitnum=87 -splitstr=""; for r in {1..87};do splitstr=$splitstr" $r"; done -for r in {1..87};do splitstr=$splitstr" $sdir/labels/seg$r/$subj.nii.gz"; done - -mirtk split-labels $num $transformed $transformedw $splitnum $splitstr - - - -#tissues -tissues="csf gm wm outlier hwm lwm" -for str in ${tissues};do -mkdir -p $sdir/labels/$str || exit 1 -done - -ln $sdir/labels/seg83/$subj.nii.gz $sdir/labels/csf -ln $sdir/labels/seg84/$subj.nii.gz $sdir/labels/outlier - -str=""; for i in ${corts};do str="$str-add $sdir/labels/seg$i/$subj.nii.gz "; done -str=`echo $str| sed -e 's:^\-add ::g'` -run mirtk calculate $str -out $sdir/labels/gm/$subj.nii.gz - -str=""; for i in ${wmcorts};do str="$str-add $sdir/labels/seg$i/$subj.nii.gz "; done -str=`echo $str| sed -e 's:^\-add ::g'` -run mirtk calculate $str -out $sdir/labels/wm/$subj.nii.gz - - -if [ -f $sdir/tissue-posteriors/hwm/$subj.nii.gz -a -f $sdir/tissue-posteriors/lwm/$subj.nii.gz ];then -run mirtk calculate $sdir/tissue-posteriors/wm/$subj.nii.gz -add $sdir/tissue-posteriors/hwm/$subj.nii.gz -add $sdir/tissue-posteriors/lwm/$subj.nii.gz -out $sdir/tissue-posteriors/$subj-wmall.nii.gz -# hwm = wm_a * hwm_t/wm_t -run mirtk calculate $sdir/tissue-posteriors/hwm/$subj.nii.gz -div-with-zero $sdir/tissue-posteriors/$subj-wmall.nii.gz -mul $sdir/labels/wm/$subj.nii.gz -out $sdir/labels/hwm/$subj.nii.gz -# lwm = wm_a * lwm_t/wm_t -run mirtk calculate $sdir/tissue-posteriors/lwm/$subj.nii.gz -div-with-zero $sdir/tissue-posteriors/$subj-wmall.nii.gz -mul $sdir/labels/wm/$subj.nii.gz -out $sdir/labels/lwm/$subj.nii.gz -# wm = wm_a - lwm - hwm -run mirtk calculate $sdir/labels/wm/$subj.nii.gz -sub $sdir/labels/hwm/$subj.nii.gz -sub $sdir/labels/lwm/$subj.nii.gz -out $sdir/labels/wm/$subj.nii.gz -rm -f $sdir/tissue-posteriors/$subj-wmall.nii.gz -fi - - - -#create MAD -if [ ! -f $sdir/MADs/$subj.nii.gz ];then -nn=5 -run mirtk calculate-gradients N4/$subj.nii.gz $sdir/MADs/$subj-grad.nii.gz 0 -run mirtk calculate-filtering $sdir/MADs/$subj-grad.nii.gz -kernel $nn -median $sdir/MADs/$subj-cur.nii.gz -run mirtk calculate $sdir/MADs/$subj-grad.nii.gz -sub $sdir/MADs/$subj-cur.nii.gz -abs -out $sdir/MADs/$subj-cur.nii.gz -run mirtk calculate-filtering $sdir/MADs/$subj-cur.nii.gz -kernel $nn -median $sdir/MADs/$subj-cur.nii.gz -run mirtk calculate $sdir/MADs/$subj-grad.nii.gz -div-with-zero $sdir/MADs/$subj-cur.nii.gz -div-with-zero 1.4826 -sq -mul 0.5 -add 1 -log -out $sdir/MADs/$subj-cur.nii.gz -run mirtk calculate N4/$subj.nii.gz -div-with-zero N4/$subj.nii.gz -mul $sdir/MADs/$subj-cur.nii.gz -add 1 -out $sdir/MADs/$subj-cur.nii.gz -run mirtk calculate $sdir/MADs/$subj-cur.nii.gz -mul 0 -add 1 -div-with-zero $sdir/MADs/$subj-cur.nii.gz -out $sdir/MADs/$subj.nii.gz -rm -f $sdir/MADs/$subj-cur.nii.gz -fi - -#create posterior penalty -str=""; for i in ${subcorts};do str="$str-add $sdir/labels/seg$i/$subj.nii.gz "; done -str=`echo $str| sed -e 's:^\-add ::g'` -run mirtk calculate $str -div-with-zero 100 -mul $sdir/MADs/$subj.nii.gz -out $sdir/MADs/$subj-subspace.nii.gz - -fi diff --git a/scripts/v1/postprocess-cortical.sh b/scripts/v1/postprocess-cortical.sh deleted file mode 100755 index fe299b8..0000000 --- a/scripts/v1/postprocess-cortical.sh +++ /dev/null @@ -1,61 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") " 1>&2; exit 1; } - -subj=$1 -sdir=segmentations-data - -run(){ - echo "$@" - "$@" || exit 1 -} - -corts=(`cat $DRAWEMDIR/parameters/cortical.csv`) -wmcorts=(`cat $DRAWEMDIR/parameters/cortical-wm.csv`) -numcorts=${#corts[*]} - -mkdir -p $sdir/cortical || exit 1 -for ((n=0;n<$numcorts;n++));do - r=${corts[$n]}; - mkdir -p $sdir/labels/seg$r-extended || exit 1 -done - - -#max prob of cortical structures + mrf regularization -segnum=0; labels=""; structs=""; -for ((n=0;n<$numcorts;n++));do - let segnum++ - r=${corts[$n]}; - w=${wmcorts[$n]}; - #merge wm,gm probs - we'll split them later using the segmentation posteriors - run mirtk calculate $sdir/labels/seg$w/$subj.nii.gz -add $sdir/labels/seg$r/$subj.nii.gz -out $sdir/labels/seg$r-extended/$subj.nii.gz; - labels="$labels $sdir/labels/seg$r-extended/$subj.nii.gz "; - structs="$structs $sdir/labels/seg$r-extended/$subj.nii.gz "; - segnumbers="$segnumbers 1 $segnum $r" -done - -run mirtk ems-hard-segmentation $segnum $labels $sdir/cortical/$subj.nii.gz -mrftimes 1 -posteriors $structs || exit 1 -run mirtk padding $sdir/cortical/$subj.nii.gz $sdir/cortical/$subj.nii.gz $sdir/cortical/$subj.nii.gz $segnumbers || exit 1 - - - - diff --git a/scripts/v1/preprocess.sh b/scripts/v1/preprocess.sh deleted file mode 100755 index b7aa50e..0000000 --- a/scripts/v1/preprocess.sh +++ /dev/null @@ -1,62 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 2 ] || { echo "usage: $(basename "$0") " 1>&2; exit 1; } -subj=$1 -age=$2 - - -if [ -n "$FSLDIR" ]; then - [ -f "$FSLDIR/bin/bet" ] || { echo "FSLDIR environment variable invalid!" 1>&2; exit 1; } -else - FSLDIR="$(cd "$(dirname "$(which bet)")"/.. && pwd)" - [ -f "$FSLDIR/bin/bet" ] || { echo "FSLDIR environment variable not set!" 1>&2; exit 1; } - export PATH="$FSLDIR/bin:$PATH" -fi - -run(){ - echo "$@" - "$@" || exit 1 -} - - -sdir=segmentations-data -dof=dofs/$subj-template-$age-n.dof.gz - -mkdir -p $sdir/brain N4 dofs bias || exit 1 - -if [ ! -f N4/$subj.nii.gz ];then - #convert image and rescale - run mirtk convert-image T2/$subj.nii.gz $sdir/brain/$subj.nii.gz -rescale 0 1000 -double - - if [ ! -f $sdir/brain/${subj}_brain_mask.nii.gz ];then - #brain extract - run bet $sdir/brain/$subj.nii.gz $sdir/brain/${subj}_brain.nii.gz -R -f 0.1 -m - fi - - #bias correct - run $DRAWEMDIR/ThirdParty/ITK/N4 3 -i $sdir/brain/$subj.nii.gz -x $sdir/brain/${subj}_brain_mask.nii.gz -o "[N4/$subj.nii.gz,bias/$subj.nii.gz]" -c "[50x50x50,0.001]" -s 2 -b "[100,3]" -t "[0.15,0.01,200]" - run mirtk calculate N4/$subj.nii.gz -mul $sdir/brain/${subj}_brain_mask.nii.gz -out N4/$subj.nii.gz - - #rescale image - run mirtk convert-image N4/$subj.nii.gz N4/$subj.nii.gz -rescale 0 1000 -double -fi diff --git a/scripts/v1/register-multi-atlas-using-gm-posteriors.sh b/scripts/v1/register-multi-atlas-using-gm-posteriors.sh deleted file mode 100755 index e6838cf..0000000 --- a/scripts/v1/register-multi-atlas-using-gm-posteriors.sh +++ /dev/null @@ -1,49 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -ge 2 -a $# -le 3 ] || { echo "usage: $(basename "$0") <#jobs>" 1>&2; exit 1; } -subj=$1 -age=$2 -njobs=1 -if [ $# -gt 2 ];then njobs=$3;fi - -run(){ - echo "$@" - "$@" || exit 1 -} - -scriptdir=$(dirname "$BASH_SOURCE") -sdir=segmentations-data -mkdir -p $sdir/gm-posteriors || exit 1 -run mirtk calculate $sdir/tissue-posteriors/gm/$subj.nii.gz -mul 100 -out $sdir/gm-posteriors/$subj.nii.gz - -mkdir -p dofs -for i in {01..20};do - atlas=ALBERT_$i - dof=dofs/$subj-$atlas-n.dof.gz - - if [ ! -f $dof ];then - run mirtk register N4/$subj.nii.gz $DRAWEMDIR/atlases/ALBERTs/T2/$atlas.nii.gz $sdir/gm-posteriors/$subj.nii.gz $DRAWEMDIR/atlases/ALBERTs/gm-posteriors-v3/$atlas.nii.gz -parin $DRAWEMDIR/parameters/ireg-multichannel-structural.cfg -dofout $dof -threads $njobs - fi - -done - diff --git a/scripts/v1/segmentation.sh b/scripts/v1/segmentation.sh deleted file mode 100755 index fb04f87..0000000 --- a/scripts/v1/segmentation.sh +++ /dev/null @@ -1,77 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") "; exit 1; } -subj=$1 - -run() -{ - echo "$@" - "$@" || exit 1 -} - -if [ ! -f segmentations/$subj-initial.nii.gz ];then -sdir=segmentations-data - -subcorts=`cat $DRAWEMDIR/parameters/subcortical-all.csv` -tissues="outlier csf gm wm hwm lwm" - - -mkdir -p segmentations $sdir/posteriors logs || exit 1; -for r in ${subcorts};do mkdir -p $sdir/posteriors/seg$r || exit 1; done -for str in ${tissues};do mkdir -p $sdir/posteriors/$str || exit 1; done - - -structs=""; saveposts=""; posts=""; num=0; -# subcortical -for r in ${subcorts};do -structs="$structs $sdir/labels/seg$r/$subj.nii.gz"; -post=$sdir/posteriors/seg$r/$subj.nii.gz -posts="$posts $post " -saveposts="$saveposts -saveprob $num $post "; -num=$(($num+1)); -done -# tissues -for str in ${tissues};do -structs="$structs $sdir/labels/$str/$subj.nii.gz"; -post=$sdir/posteriors/$str/$subj.nii.gz -posts="$posts $post " -saveposts="$saveposts -saveprob $num $post "; -num=$(($num+1)); -done - - - -# segmentation -run mirtk draw-em N4/$subj.nii.gz 27 $structs segmentations/$subj-em.nii.gz -padding 0 -mrf $DRAWEMDIR/parameters/connectivities.mrf -tissues 1 21 1 22 1 23 3 24 25 26 -hui -postpenalty $sdir/MADs/$subj-subspace.nii.gz $saveposts 1>logs/$subj-em 2>logs/$subj-em-err - -# add hwm and lwm posterior probability to wm -run mirtk calculate $sdir/posteriors/hwm/$subj.nii.gz -add $sdir/posteriors/lwm/$subj.nii.gz -add $sdir/posteriors/wm/$subj.nii.gz -out $sdir/posteriors/wm/$subj.nii.gz - -# posteriors [0,1] -> [0,100] -for post in ${posts};do -run mirtk calculate $post -mul 100 -out $post -done - -# set whole wm to 2000 and whole gm to 1000 and final label numbers for the rest (subcortical, csf, out) -run mirtk padding segmentations/$subj-em.nii.gz segmentations/$subj-em.nii.gz segmentations/$subj-initial.nii.gz $DRAWEMDIR/parameters/seg-numbers.csv -fi diff --git a/scripts/v1/separate-hemispheres.sh b/scripts/v1/separate-hemispheres.sh deleted file mode 100755 index 5c90c59..0000000 --- a/scripts/v1/separate-hemispheres.sh +++ /dev/null @@ -1,94 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -eq 1 ] || { echo "usage: $(basename "$0") "; exit 1; } -subj=$1 - -run(){ - echo "$@" - "$@" || exit 1 -} - -scriptdir=$(dirname "$BASH_SOURCE") - - -if [ ! -f segmentations/${subj}_L_white.nii.gz -o ! -f segmentations/${subj}_R_white.nii.gz -o ! -f segmentations/${subj}_L_pial.nii.gz -o ! -f segmentations/${subj}_R_pial.nii.gz ];then - -# create initial labels files -suffix=-sephemi -$scriptdir/postprocess.sh $subj $suffix - - -# structures for the left, right hemisphere and those that extend to both -left=`cat $DRAWEMDIR/parameters/structures-left.csv` -right=`cat $DRAWEMDIR/parameters/structures-right.csv` -both=`cat $DRAWEMDIR/parameters/structures-both.csv` -numleft=`echo $left|wc -w` -numright=`echo $right|wc -w` -numboth=`echo $both|wc -w` - -# left, right, both hemispheres -run mirtk calculate segmentations/"$subj"_labels$suffix.nii.gz -mul 0 -out segmentations/$subj-empty.nii.gz -run mirtk padding segmentations/$subj-empty.nii.gz segmentations/"$subj"_labels$suffix.nii.gz segmentations/$subj-L-hemisphere.nii.gz $numleft `echo $left` 1 -run mirtk padding segmentations/$subj-empty.nii.gz segmentations/"$subj"_labels$suffix.nii.gz segmentations/$subj-R-hemisphere.nii.gz $numright `echo $right` 1 -run mirtk padding segmentations/$subj-empty.nii.gz segmentations/"$subj"_labels$suffix.nii.gz segmentations/$subj-both-hemisphere.nii.gz $numboth `echo $both` 1 -run mirtk padding segmentations/$subj-both-hemisphere.nii.gz segmentations/"$subj"_all_labels$suffix.nii.gz segmentations/$subj-both-hemisphere.nii.gz 85 1 - -# compute a cutting plane based on dmap -run mirtk calculate-distance-map segmentations/$subj-L-hemisphere.nii.gz segmentations/$subj-L-hemisphere-dmap.nii.gz -run mirtk calculate-distance-map segmentations/$subj-R-hemisphere.nii.gz segmentations/$subj-R-hemisphere-dmap.nii.gz -run mirtk smooth-image segmentations/$subj-L-hemisphere-dmap.nii.gz segmentations/$subj-L-hemisphere-dmap.nii.gz 1 -float -run mirtk smooth-image segmentations/$subj-R-hemisphere-dmap.nii.gz segmentations/$subj-R-hemisphere-dmap.nii.gz 1 -float -run mirtk calculate segmentations/$subj-R-hemisphere-dmap.nii.gz -sub segmentations/$subj-L-hemisphere-dmap.nii.gz -mask-below 0 -inside 1 -outside 0 -out segmentations/$subj-hemisphere-cut.nii.gz - -# split structures that are in both hemispheres into left and right according to the cutting plane -run mirtk calculate segmentations/$subj-hemisphere-cut.nii.gz -mul segmentations/$subj-both-hemisphere.nii.gz -add segmentations/$subj-L-hemisphere.nii.gz -out segmentations/$subj-L-hemisphere.nii.gz -run mirtk calculate segmentations/$subj-hemisphere-cut.nii.gz -sub 1 -mul -1 -mul segmentations/$subj-both-hemisphere.nii.gz -add segmentations/$subj-R-hemisphere.nii.gz -out segmentations/$subj-R-hemisphere.nii.gz -rm segmentations/$subj-both-hemisphere.nii.gz segmentations/$subj-hemisphere-cut.nii.gz segmentations/$subj-L-hemisphere-dmap.nii.gz segmentations/$subj-R-hemisphere-dmap.nii.gz segmentations/$subj-empty.nii.gz - -rmlabelswhite="5 1 2 4 6 8" -rmlabelspial="4 1 4 6 8" - -for h in L R;do - for surf in white pial;do - # left, right wm+dgm / gm+wm+cgm - rmlabels=rmlabels$surf - run mirtk padding segmentations/$subj-$h-hemisphere.nii.gz segmentations/"$subj"_tissue_labels$suffix.nii.gz segmentations/${subj}_${h}_${surf}_init.nii.gz ${!rmlabels} 0 - - # keep only connected components with volume > 5% volume of first component - $scriptdir/clear-small-components.sh segmentations/${subj}_${h}_${surf}_init.nii.gz segmentations/${subj}_${h}_${surf}_unfilled.nii.gz - - # fill holes - run mirtk fill-holes segmentations/${subj}_${h}_${surf}_unfilled.nii.gz segmentations/${subj}_${h}_${surf}.nii.gz - - # clean up - rm segmentations/${subj}_${h}_${surf}_init.nii.gz segmentations/${subj}_${h}_${surf}_unfilled.nii.gz - done - # clean up - rm segmentations/$subj-$h-hemisphere.nii.gz -done - -# clean up -rm segmentations/"$subj"_all_labels$suffix.nii.gz segmentations/"$subj"_labels$suffix.nii.gz segmentations/"$subj"_tissue_labels$suffix.nii.gz - -fi - - diff --git a/scripts/v1/tissue-priors.sh b/scripts/v1/tissue-priors.sh deleted file mode 100755 index 9608cb7..0000000 --- a/scripts/v1/tissue-priors.sh +++ /dev/null @@ -1,72 +0,0 @@ -#!/bin/bash -# ============================================================================ -# Developing brain Region Annotation With Expectation-Maximization (Draw-EM) -# -# Copyright 2013-2016 Imperial College London -# Copyright 2013-2016 Andreas Schuh -# Copyright 2013-2016 Antonios Makropoulos -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ============================================================================ - - -[ $# -ge 3 ] || { echo "usage: $(basename "$0") [<#jobs>]" 1>&2; exit 1; } -subj=$1 -age=$2 -atlasname=$3 -njobs=1 -if [ $# -gt 3 ];then njobs=$4;fi - -run(){ - echo "$@" - "$@" || exit 1 -} - - -#registration -dof=dofs/$subj-template-$age-n.dof.gz -mkdir -p $(dirname "${dof}") - -if [ ! -f $dof ];then - run mirtk register N4/$subj.nii.gz $DRAWEMDIR/atlases/$atlasname/T2/template-$age.nii.gz -dofout $dof -parin $DRAWEMDIR/parameters/ireg.cfg -threads $njobs -fi - -sdir=segmentations-data -if [ ! -f $sdir/tissue-initial-segmentations/$subj.nii.gz ];then - -echo "creating $subj tissue priors" - -mkdir -p $sdir $sdir/template $sdir/tissue-posteriors $sdir/tissue-initial-segmentations || exit 1 -structures="csf gm wm outlier ventricles cerebstem dgm hwm lwm" -for str in ${structures};do -mkdir -p $sdir/template/$str $sdir/tissue-posteriors/$str || exit 1 -done - - -strnum=0 -emsstructures="" -emsposts="" - -for str in ${structures};do -emsposts="$emsposts -saveprob $strnum $sdir/tissue-posteriors/$str/$subj.nii.gz" -strnum=$(($strnum+1)) - -strems=$sdir/template/$str/$subj.nii.gz -run mirtk transform-image $DRAWEMDIR/atlases/$atlasname/atlas-9/structure$strnum/$age.nii.gz $strems -dofin $dof -target N4/$subj.nii.gz -interp Linear -emsstructures="$emsstructures $strems" -done - -mkdir -p logs -run mirtk draw-em N4/$subj.nii.gz 9 $emsstructures $sdir/tissue-initial-segmentations/$subj.nii.gz -padding 0 -mrf $DRAWEMDIR/parameters/conn_tissues_ven_cstem_dgm_hwm_lwm.mrf -tissues 1 3 1 0 1 1 3 2 7 8 -hui -relaxtimes 2 $emsposts 1>logs/$subj-tissue-em 2>logs/$subj-tissue-em-err - -fi diff --git a/segmentation.png b/segmentation.png new file mode 100644 index 0000000..d69eb41 Binary files /dev/null and b/segmentation.png differ diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 58be792..f0cef2d 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -47,3 +47,5 @@ add_drawem_command(em-hard-segmentation) add_drawem_command(kmeans) add_drawem_command(normalize) add_drawem_command(split-labels) + +mirtk_add_executable(neonatal-segmentation) diff --git a/tools/draw-em.cc b/tools/draw-em.cc index 41c2641..5641281 100644 --- a/tools/draw-em.cc +++ b/tools/draw-em.cc @@ -576,9 +576,9 @@ int main(int argc, char **argv) clock_t end = clock(); - int elapsed_secs = iround( double(end - begin) / CLOCKS_PER_SEC); - int elapsed_mins = elapsed_secs/60; - elapsed_secs%=elapsed_mins; + int elapsed_secs = round( double(end - begin) / CLOCKS_PER_SEC); + int elapsed_mins = round(elapsed_secs/60); + elapsed_secs = round(elapsed_secs%60); std::cout<<"elapsed time: "< The directory used to run the script and output the files. -c / -cleanup <0/1> Whether cleanup of temporary files is required (default: 1) -p / -save-posteriors <0/1> Whether the structures' posteriors are required (default: 0) + -atlas Atlas used for the tissue priors (default: non-rigid-v2) -t / -threads Number of threads (CPU cores) allowed for the registration to run in parallel (default: 1) -v / -verbose <0/1> Whether the script progress is reported (default: 1) -h / -help / --help Print usage. @@ -65,15 +66,16 @@ datadir=`pwd` posteriors=0 # whether to output posterior probability maps threads=1 verbose=1 +command="$@" atlasname=non-rigid-v2 -version="v1.1" while [ $# -gt 0 ]; do case "$3" in -c|-cleanup) shift; cleanup=$3; ;; -d|-data-dir) shift; datadir=$3; ;; -p|-save-posteriors) shift; posteriors=$3; ;; + -atlas) shift; atlasname=$3; ;; -t|-threads) shift; threads=$3; ;; -v|-verbose) shift; verbose=$3; ;; -h|-help|--help) usage; ;; @@ -91,8 +93,10 @@ else fi cd $datadir +version=`git -C "$DRAWEMDIR" branch | grep \* | cut -d ' ' -f2` +gitversion=`git -C "$DRAWEMDIR" rev-parse HEAD` -[ $verbose -le 0 ] || { echo "DrawEM multi atlas $version +[ $verbose -le 0 ] || { echo "DrawEM multi atlas $version (branch version: $gitversion) Subject: $subj Age: $age Directory: $datadir @@ -100,43 +104,55 @@ Posteriors: $posteriors Cleanup: $cleanup Threads: $threads -$BASH_SOURCE $@ +$BASH_SOURCE $command ----------------------------"; } mkdir -p logs || exit 1 +# log function run() { - [ $verbose -le 0 ] || echo -n "$@..." - "$DRAWEMDIR/scripts/$version/$@" 1>>logs/$subj 2>>logs/$subj-err - if [ $? -eq 0 ]; then - [ $verbose -le 0 ] || echo " done" - else - [ $verbose -le 0 ] || echo " failed: see log file logs/$subj-err for details" + echo "$@" + "$@" + if [ ! $? -eq 0 ]; then + echo "$@ : failed" + exit 1 + fi +} + +# make run function global +typeset -fx run + +run_script() +{ + echo "$@" + "$DRAWEMDIR/scripts/$@" + if [ ! $? -eq 0 ]; then + echo "$DRAWEMDIR/scripts/$@ : failed" exit 1 fi } rm -f logs/$subj logs/$subj-err -run preprocess.sh $subj +run_script preprocess.sh $subj # phase 1 tissue segmentation -run tissue-priors.sh $subj $age $atlasname $threads +run_script tissue-priors.sh $subj $age $atlasname $threads # registration using gm posterior + image -run register-multi-atlas-using-gm-posteriors.sh $subj $age $threads +run_script register-multi-atlas-using-gm-posteriors.sh $subj $age $threads # structural segmentation -run labels-multi-atlas.sh $subj -run segmentation.sh $subj +run_script labels-multi-atlas.sh $subj +run_script segmentation.sh $subj # post-processing -run separate-hemispheres.sh $subj -run correct-segmentation.sh $subj -run postprocess.sh $subj +run_script separate-hemispheres.sh $subj +run_script correct-segmentation.sh $subj +run_script postprocess.sh $subj # if probability maps are required -[ "$posteriors" == "0" -o "$posteriors" == "no" -o "$posteriors" == "false" ] || run postprocess-pmaps.sh $subj +[ "$posteriors" == "0" -o "$posteriors" == "no" -o "$posteriors" == "false" ] || run_script postprocess-pmaps.sh $subj # cleanup if [ "$cleanup" == "1" -o "$cleanup" == "yes" -o "$cleanup" == "true" ] && [ -f "segmentations/${subj}_labels.nii.gz" ];then - run clear-data.sh $subj + run_script clear-data.sh $subj rm -f logs/$subj logs/$subj-err rmdir logs 2> /dev/null # may fail if other log files exist fi