Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

ENH: Added BRAINSABC,BCD,BRAINSFit,BRAINSCommon

  • Loading branch information...
commit d98591cbcc94919fea3490265fa41ff3efc0f3e2 1 parent 339022c
@hjmjohnson hjmjohnson authored
Showing with 42,293 additions and 0 deletions.
  1. +41 −0 BRAINSABC/CMakeLists.txt
  2. +14 −0 BRAINSABC/CTestConfig.cmake
  3. +1 −0  BRAINSABC/TestSuite/AtlasPVDefinition.xml.md5
  4. +23 −0 BRAINSABC/TestSuite/BRAINSABCTest.cxx
  5. +54 −0 BRAINSABC/TestSuite/BlendImageFilterTest.cxx
  6. +60 −0 BRAINSABC/TestSuite/CMakeLists.txt
  7. +347 −0 BRAINSABC/TestSuite/simpleEM.cxx
  8. +143 −0 BRAINSABC/TestSuite/simpleEM.xml
  9. +1 −0  BRAINSABC/TestSuite/test_data/000517510_000517510_lmks.fcsv.md5
  10. +1 −0  BRAINSABC/TestSuite/test_data/EMS-Param.xml.md5
  11. +1 −0  BRAINSABC/TestSuite/test_data/EMS.log.md5
  12. +1 −0  BRAINSABC/TestSuite/test_data/EMS.xml.md5
  13. +1 −0  BRAINSABC/TestSuite/test_data/FAST_logger.log.md5
  14. +1 −0  BRAINSABC/TestSuite/test_data/ISO_T1_REP0.nii.gz.md5
  15. +1 −0  BRAINSABC/TestSuite/test_data/ISO_T1_REP0_to_template_EMS.bspline.md5
  16. +1 −0  BRAINSABC/TestSuite/test_data/ISO_T1_REP1.nii.gz.md5
  17. +1 −0  BRAINSABC/TestSuite/test_data/ISO_T2_REP0.nii.gz.md5
  18. +1 −0  BRAINSABC/TestSuite/test_data/T1_REP0.nii.gz.md5
  19. +1 −0  BRAINSABC/TestSuite/test_data/T1_REP1.nii.gz.md5
  20. +1 −0  BRAINSABC/TestSuite/test_data/T2_REP0.nii.gz.md5
  21. +1 −0  BRAINSABC/TestSuite/test_data/binary_release.test.md5
  22. +1 −0  BRAINSABC/TestSuite/test_data/medium_T1_REP0.nii.gz.md5
  23. +1 −0  BRAINSABC/TestSuite/test_data/medium_T1_REP1.nii.gz.md5
  24. +1 −0  BRAINSABC/TestSuite/test_data/medium_T2_REP0.nii.gz.md5
  25. +1 −0  BRAINSABC/TestSuite/test_data/qqq.logger.md5
  26. +1 −0  BRAINSABC/TestSuite/test_data/refResults/labels.nii.gz.md5
  27. +1 −0  BRAINSABC/TestSuite/test_data/refResults/medium_T1_REP0_labels_BRAINSABC.nii.gz.md5
  28. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_corrected_BRAINSABC.nii.gz.md5
  29. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_corrected_BRAINSABC.nrrd.md5
  30. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_labels_BRAINSABC.nii.gz.md5
  31. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_labels_BRAINSABC.nrrd.md5
  32. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_posterior0_BRAINSABC.nrrd.md5
  33. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_posterior1_BRAINSABC.nrrd.md5
  34. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_posterior2_BRAINSABC.nrrd.md5
  35. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_registered_BRAINSABC.nrrd.md5
  36. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_template_affine_BRAINSABC.nrrd.md5
  37. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_template_warped_BRAINSABC.nrrd.md5
  38. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP1_corrected_BRAINSABC.nrrd.md5
  39. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP1_registered_BRAINSABC.nrrd.md5
  40. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T2_REP0_corrected_BRAINSABC.nrrd.md5
  41. +1 −0  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T2_REP0_registered_BRAINSABC.nrrd.md5
  42. +1 −0  BRAINSABC/TestSuite/test_data/run_test.sh.md5
  43. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0.nii.gz.md5
  44. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_small_ISO_T1_REP1_EMS.affine.md5
  45. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_small_ISO_T2_REP0_EMS.affine.md5
  46. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_template_EMS.affine.md5
  47. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_template_EMS.bspline.md5
  48. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP1.nii.gz.md5
  49. +1 −0  BRAINSABC/TestSuite/test_data/small_ISO_T2_REP0.nii.gz.md5
  50. +1 −0  BRAINSABC/TestSuite/test_data/small_T1_REP0.nii.gz.md5
  51. +1 −0  BRAINSABC/TestSuite/test_data/small_T1_REP0_to_template_EMS.bspline.md5
  52. +1 −0  BRAINSABC/TestSuite/test_data/small_T1_REP1.nii.gz.md5
  53. +1 −0  BRAINSABC/TestSuite/test_data/small_T2_REP0.nii.gz.md5
  54. +1 −0  BRAINSABC/TestSuite/test_data/test.sh.md5
  55. +1 −0  BRAINSABC/TestSuite/test_data/tmp_load.b2.md5
  56. +122 −0 BRAINSABC/brainseg/AtlasCropImageSource.h
  57. +291 −0 BRAINSABC/brainseg/AtlasCropImageSource.txx
  58. +384 −0 BRAINSABC/brainseg/AtlasDefinition.cxx
  59. +361 −0 BRAINSABC/brainseg/AtlasDefinition.h
  60. +195 −0 BRAINSABC/brainseg/AtlasRegistrationMethod.h
  61. +772 −0 BRAINSABC/brainseg/AtlasRegistrationMethod.txx
  62. +3 −0  BRAINSABC/brainseg/AtlasRegistrationMethod_float+float.cxx
  63. +1,382 −0 BRAINSABC/brainseg/BRAINSABC.cxx
  64. +271 −0 BRAINSABC/brainseg/BRAINSABC.xml
  65. +98 −0 BRAINSABC/brainseg/BRAINSABCUtilities.cxx
  66. +106 −0 BRAINSABC/brainseg/BRAINSABCUtilities.h
  67. +294 −0 BRAINSABC/brainseg/BRAINSABCUtilities.txx
  68. +15 −0 BRAINSABC/brainseg/BRAINSABCimple.xml
  69. +73 −0 BRAINSABC/brainseg/BRAINSCleanMask.cxx
  70. +39 −0 BRAINSABC/brainseg/BRAINSCleanMask.xml
  71. +103 −0 BRAINSABC/brainseg/CMakeLists.txt
  72. +231 −0 BRAINSABC/brainseg/ComputeDistributions.h
  73. +127 −0 BRAINSABC/brainseg/DenoiseFiltering.h
  74. +154 −0 BRAINSABC/brainseg/EMSParameters.cxx
  75. +132 −0 BRAINSABC/brainseg/EMSParameters.h
  76. +324 −0 BRAINSABC/brainseg/EMSegmentationFilter.h
  77. +1,839 −0 BRAINSABC/brainseg/EMSegmentationFilter.txx
  78. +7 −0 BRAINSABC/brainseg/EMSegmentationFilter_float+float.cxx
  79. +393 −0 BRAINSABC/brainseg/EMSplitPriorMST.h
  80. +58 −0 BRAINSABC/brainseg/ESLR.cxx
  81. +80 −0 BRAINSABC/brainseg/ESLR.xml
  82. +132 −0 BRAINSABC/brainseg/ExtractSingleLargestRegion.cxx
  83. +23 −0 BRAINSABC/brainseg/ExtractSingleLargestRegion.h
  84. +133 −0 BRAINSABC/brainseg/GenerateLabelMapFromProbabilityMap.cxx
  85. +47 −0 BRAINSABC/brainseg/GenerateLabelMapFromProbabilityMap.xml
  86. +185 −0 BRAINSABC/brainseg/LLSBiasCorrector.h
  87. +867 −0 BRAINSABC/brainseg/LLSBiasCorrector.txx
  88. +41 −0 BRAINSABC/brainseg/MSTEdge.h
  89. +443 −0 BRAINSABC/brainseg/QHullMSTClusteringProcess.cxx
  90. +73 −0 BRAINSABC/brainseg/QHullMSTClusteringProcess.h
  91. +64 −0 BRAINSABC/brainseg/RegistrationParameters.h
  92. +58 −0 BRAINSABC/brainseg/StandardizeMaskIntensity.cxx
  93. +134 −0 BRAINSABC/brainseg/StandardizeMaskIntensity.h
  94. +55 −0 BRAINSABC/brainseg/filterFloatImages.cxx
  95. +13 −0 BRAINSABC/brainseg/filterFloatImages.h
  96. +111 −0 BRAINSABC/brainseg/itkBlendImageFilter.h
  97. +88 −0 BRAINSABC/brainseg/itkBlendImageFilter.txx
  98. +24 −0 BRAINSABC/brainseg/normalA-30.xml
  99. +236 −0 BRAINSABC/brainseg/oldCode.cxx
  100. +81 −0 BRAINSABC/common/Heap.h
  101. +265 −0 BRAINSABC/common/Heap.txx
  102. +102 −0 BRAINSABC/common/Log.cxx
  103. +76 −0 BRAINSABC/common/Log.h
  104. +17 −0 BRAINSABC/common/mu.h
  105. +47 −0 BRAINSABC/common/muCygwin.h
  106. +68 −0 BRAINSABC/common/muException.h
  107. +46 −0 BRAINSABC/common/muMacros.h
  108. +66 −0 BRAINSABC/common/muUtils.h
  109. +17 −0 BRAINSABC/common/testexc.cxx
  110. +89 −0 BRAINSABC/common/testheap.cxx
  111. +26 −0 BRAINSABC/common/testmacro.cxx
  112. +85 −0 BRAINSABC/notes.BRAINSABC
  113. +43 −0 BRAINSABC/qhull/CMakeLists.txt
  114. +1,085 −0 BRAINSABC/qhull/Changes.txt
  115. +97 −0 BRAINSABC/qhull/Makefile.am
  116. +697 −0 BRAINSABC/qhull/Makefile.in
  117. +191 −0 BRAINSABC/qhull/Makefile.txt
  118. +1,507 −0 BRAINSABC/qhull/geom.c
  119. +225 −0 BRAINSABC/qhull/geom.h
  120. +2,600 −0 BRAINSABC/qhull/geom2.c
  121. +2,591 −0 BRAINSABC/qhull/global.c
  122. +5,346 −0 BRAINSABC/qhull/io.c
  123. +212 −0 BRAINSABC/qhull/io.h
  124. +535 −0 BRAINSABC/qhull/mem.c
  125. +352 −0 BRAINSABC/qhull/mem.h
  126. +4,330 −0 BRAINSABC/qhull/merge.c
  127. +300 −0 BRAINSABC/qhull/merge.h
  128. +1,431 −0 BRAINSABC/qhull/poly.c
  129. +361 −0 BRAINSABC/qhull/poly.h
  130. +3,888 −0 BRAINSABC/qhull/poly2.c
  131. +364 −0 BRAINSABC/qhull/qconvex.c
  132. +355 −0 BRAINSABC/qhull/qdelaun.c
  133. +359 −0 BRAINSABC/qhull/qhalf.c
  134. +1,808 −0 BRAINSABC/qhull/qhull.c
  135. +1,165 −0 BRAINSABC/qhull/qhull.h
  136. +155 −0 BRAINSABC/qhull/qhull_a.h
Sorry, we could not display the entire diff because too many files (709) changed.
View
41 BRAINSABC/CMakeLists.txt
@@ -0,0 +1,41 @@
+project(BRAINSABC)
+set(LOCAL_PROJECT_NAME BRAINSABC)
+cmake_minimum_required(VERSION 2.8)
+cmake_policy(VERSION 2.8)
+
+enable_testing()
+include(CTest)
+
+find_package(BRAINSCommonLib NO_MODULE REQUIRED)
+include(${BRAINSCommonLib_USE_FILE})
+
+include(${BRAINSCommonLib_BUILDSCRIPTS_DIR}/PreventInSourceBuilds.cmake)
+include(${BRAINSCommonLib_BUILDSCRIPTS_DIR}/CMakeBuildMacros.cmake)
+include(${BRAINSCommonLib_BUILDSCRIPTS_DIR}/SEMMacroBuildCLI.cmake)
+include(${BRAINSCommonLib_BUILDSCRIPTS_DIR}/CMakeBRAINS3BuildMacros.cmake)
+include(${BRAINSCommonLib_BUILDSCRIPTS_DIR}/IJMacros.txt)
+
+###
+SETIFEMPTY(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/lib)
+SETIFEMPTY(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/lib)
+SETIFEMPTY(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/bin)
+SETIFEMPTY(CMAKE_BUNDLE_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/bin)
+link_directories(${CMAKE_LIBRARY_OUTPUT_DIRECTORY} ${CMAKE_ARCHIVE_OUTPUT_DIRECTORY})
+
+if(NOT ITK_FOUND)
+ find_package(ITK REQUIRED)
+ include(${ITK_USE_FILE})
+endif(NOT ITK_FOUND)
+
+#-----------------------------------------------------------------------------
+# Output directories.
+#
+#SETOPTIONALDEBUGIMAGEVIEWER()
+
+###
+add_subdirectory(brainseg)
+
+if(1)
+ add_subdirectory(TestSuite)
+endif(1)
+
View
14 BRAINSABC/CTestConfig.cmake
@@ -0,0 +1,14 @@
+## This file should be placed in the root directory of your project.
+## Then modify the CMakeLists.txt file in the root directory of your
+## project to incorporate the testing dashboard.
+## # The following are required to uses Dart and the Cdash dashboard
+## enable_testing()
+## include(CTest)
+set(CTEST_PROJECT_NAME "BRAINSABC")
+set(CTEST_NIGHTLY_START_TIME "00:00:00 EST")
+
+set(CTEST_DROP_METHOD "http")
+set(CTEST_DROP_SITE "testing.psychiatry.uiowa.edu")
+set(CTEST_DROP_LOCATION "/CDash/submit.php?project=BRAINSABC")
+set(CTEST_DROP_SITE_CDASH TRUE)
+set(CTEST_TEST_TIMEOUT 3600) ## Set timeout to one hour for now. There is a lot of work to be done.
View
1  BRAINSABC/TestSuite/AtlasPVDefinition.xml.md5
@@ -0,0 +1 @@
+24d6e8d67a2fa417853e1b86449c8cfb
View
23 BRAINSABC/TestSuite/BRAINSABCTest.cxx
@@ -0,0 +1,23 @@
+//
+//A test driver to append the
+//itk image processing test
+//commands to an
+//the SEM compatibile program
+//
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#ifdef WIN32
+#define MODULE_IMPORT __declspec(dllimport)
+#else
+#define MODULE_IMPORT
+#endif
+
+extern "C" MODULE_IMPORT int ModuleEntryPoint(int, char* []);
+
+int BRAINSABCTest(int argc, char** argv)
+{
+ return ModuleEntryPoint(argc, argv);
+}
+
View
54 BRAINSABC/TestSuite/BlendImageFilterTest.cxx
@@ -0,0 +1,54 @@
+#include "itkImage.h"
+#include "itkBlendImageFilter.h"
+#include "itkRandomImageSource.h"
+#include "itkImageRegionConstIterator.h"
+#include "vnl/vnl_math.h"
+
+int main(int, char * *)
+{
+ typedef itk::Image<float, 2> ImageType;
+
+ std::cout << "Create input image using RandomImageSource" << std::endl;
+ ImageType::Pointer images[2];
+ for( unsigned i = 0; i < 2; i++ )
+ {
+ typedef itk::RandomImageSource<ImageType> SourceType;
+ SourceType::Pointer source = SourceType::New();
+ ImageType::SizeValueType size[2] = {64, 64};
+ source->SetSize( size );
+ source->SetMin(0.0);
+ source->SetMax(1.0);
+ source->Update();
+ images[i] = source->GetOutput();
+ }
+ typedef itk::BlendImageFilter<ImageType, ImageType>
+ BlendImageFilterType;
+ BlendImageFilterType::Pointer filter =
+ BlendImageFilterType::New();
+ filter->SetInput1(images[0]);
+ filter->SetInput2(images[1]);
+ filter->SetBlend1(0.2);
+ filter->SetBlend2(0.8);
+ filter->Update();
+ ImageType::Pointer blendImage = filter->GetOutput();
+ itk::ImageRegionConstIterator<ImageType>
+ it1(images[0],
+ images[0]->GetLargestPossibleRegion() ),
+ it2(images[1],
+ images[1]->GetLargestPossibleRegion() ),
+ itBlend(blendImage,
+ blendImage->GetLargestPossibleRegion() );
+ for( ; !it1.IsAtEnd() && !it2.IsAtEnd() && !itBlend.IsAtEnd();
+ ++it1, ++it2, ++itBlend )
+ {
+ float blend = (it1.Get() * 0.2) + (it2.Get() * 0.8);
+ if( vcl_fabs(blend - itBlend.Get() ) > 0.0001 )
+ {
+ std::cerr << "Expected " << blend << " found " << itBlend.Get()
+ << std::endl;
+ exit(1);
+ }
+ }
+ exit(0);
+}
+
View
60 BRAINSABC/TestSuite/CMakeLists.txt
@@ -0,0 +1,60 @@
+list(APPEND ExternalData_URL_TEMPLATES
+ # Local data store populated by the ITK pre-commit hook
+ "file:///${CMAKE_SOURCE_DIR}/.ExternalData/%(algo)/%(hash)"
+ # Data published by Iowa Psychiatry web interface
+ "http://www.psychiatry.uiowa.edu/users/brainstestdata/ctestdata/%(algo)/%(hash)"
+
+ # Data published by MIDAS
+ "http://midas.kitware.com/api/rest/midas.bitstream.by.hash?hash=%(hash)&algorithm=%(algo)"
+
+ # Data published by developers using git-gerrit-push.
+ "http://www.itk.org/files/ExternalData/%(algo)/%(hash)"
+ )
+
+# Tell ExternalData commands to transform raw files to content links.
+# TODO: Condition this feature on presence of our pre-commit hook.
+set(ExternalData_LINK_CONTENT MD5)
+include(${BRAINSCommonLib_BUILDSCRIPTS_DIR}/ExternalData.cmake)
+
+include_directories(
+ ${BRAINSABC_SOURCE_DIR}/brainseg
+ ${BRAINSABC_SOURCE_DIR}/common
+ ${BRAINSABC_BINARY_DIR}/brainseg
+)
+
+
+
+MakeTestDriverFromSEMTool(BRAINSABC BRAINSABCTest.cxx)
+
+ExternalData_add_test( FetchData NAME BRAINSABCLongTest
+ COMMAND ${LAUNCH_EXE} $<TARGET_FILE:BRAINSABCTestDriver>
+ --compare DATA{test_data/refResults/labels.nii.gz}
+ ${CMAKE_CURRENT_BINARY_DIR}/labels.nii.gz
+ --compareIntensityTolerance 1
+ --compareRadiusTolerance 1
+ --compareNumberOfPixelsTolerance 10000
+ BRAINSABCTest
+ --inputVolumes DATA{test_data/small_ISO_T1_REP0.nii.gz} --inputVolumes DATA{test_data/small_ISO_T2_REP0.nii.gz}
+ --outputVolumes ${CMAKE_CURRENT_BINARY_DIR}/T1_cleaned.nii.gz,${CMAKE_CURRENT_BINARY_DIR}/T2_cleaned.nii.gz
+ --outputLabels ${CMAKE_CURRENT_BINARY_DIR}/labels.nii.gz
+ --outputDirtyLabels ${CMAKE_CURRENT_BINARY_DIR}/dirty_labels.nii.gz
+ --posteriorTemplate ${CMAKE_CURRENT_BINARY_DIR}/POST_%s.nii.gz
+ --inputVolumeTypes T1,T2
+ --filterIteration 3
+ --maxIterations 2
+ --maxBiasDegree 2
+ --debuglevel 0
+ --outputFormat NIFTI
+ --outputDir ${CMAKE_CURRENT_BINARY_DIR}
+ --gridSize 5,3,4
+ --atlasDefinition ${CMAKE_CURRENT_BINARY_DIR}/../../src/bin/Atlas/Atlas_20110607/AtlasPVDefinition.xml
+)
+set_tests_properties(BRAINSABCLongTest PROPERTIES TIMEOUT 6500)
+
+if(0) # This should be restored after fixing.
+add_executable(BlendImageFilterTest BlendImageFilterTest.cxx)
+target_link_libraries(BlendImageFilterTest ${ITK_LIBRARIES})
+ExternalData_add_test( FetchData NAME BlendImageFilterTest COMMAND ${LAUNCH_EXE} $<TARGET_FILE:BlendImageFilterTest> )
+endif()
+
+ExternalData_Add_Target( FetchData ) # Name of data management target
View
347 BRAINSABC/TestSuite/simpleEM.cxx
@@ -0,0 +1,347 @@
+#include "simpleEMCLP.h"
+
+#include <itksys/SystemTools.hxx>
+
+#include "itkCastImageFilter.h"
+#include "itkImage.h"
+#include "itkImageFileWriter.h"
+#include "itkImageFileReader.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "itkNumericTraits.h"
+#include "itkMultiplyImageFilter.h"
+#include "itkResampleImageFilter.h"
+#include "itkRescaleIntensityImageFilter.h"
+#include "itkVersion.h"
+#include "AtlasCropImageSource.h"
+
+#include <vector>
+
+// Use manually instantiated classes for the big program chunks
+// #define MU_MANUAL_INSTANTIATION
+#include "EMSegmentationFilter.h"
+#include "itkLargestForegroundFilledMaskImageFilter.h"
+
+// #undef MU_MANUAL_INSTANTIATION
+
+#include "filterFloatImages.h"
+
+#include <iostream>
+#include <string>
+#include <sstream>
+
+#include <stdlib.h>
+// template <class inputPixelType,class outputPixelType,class priorPixelType>
+template <class priorPixelType, class inputPixelType>
+int simpleRunEMS( std::string t1Volume,
+ std::string t2Volume,
+ std::string pdVolume,
+ std::string templateVolume,
+ int maxIteration,
+ std::string OutputFileNamePrefix,
+ std::vector<std::string> & priorsList,
+ std::vector<double> & priorsWeightList,
+ bool warp,
+ int degreeOfBiasFieldCorrection,
+ double likelihoodTolerance)
+{
+ int status = -1;
+ // PARSE_ARGS;
+ /** read parameters **/
+ /** - read target input images **/
+ std::string T1FileName(t1Volume);
+ std::string T2FileName(t2Volume);
+ std::string PDFileName(pdVolume);
+
+ typedef itk::Image<float, 3> FloatImageType;
+ typedef itk::Image<unsigned char, 3> ByteImageType;
+ typedef itk::Image<short, 3> ShortImageType;
+
+ typedef FloatImageType::Pointer FloatImagePointer;
+ typedef ByteImageType::Pointer ByteImagePointer;
+ typedef ShortImageType::Pointer ShortImagePointer;
+
+ typedef itk::Image<inputPixelType, 3> InputImageType;
+ typedef itk::ImageFileReader<InputImageType> InputImageReaderType;
+ typedef itk::RescaleIntensityImageFilter<InputImageType, InputImageType> RescaleImageFilterType;
+ typedef typename InputImageType::Pointer InputImagePointer;
+
+ typedef itk::ImageFileReader<InputImageType> templateReaderType;
+ typename templateReaderType::Pointer templateReader = templateReaderType::New();
+ templateReader->SetFileName(templateVolume);
+ templateReader->Update();
+ std::vector<InputImagePointer> images;
+
+ /* IMAGES MUST BE Non-Negative or the results will fail. */
+ std::vector<std::string> inputImageFilenames;
+ if( T1FileName == "" )
+ {
+ std::cerr << " T1 file is missing, other modality will be used as a reference image.\n ";
+ }
+ else
+ {
+ typename InputImageReaderType::Pointer inputReaderT1 = InputImageReaderType::New();
+ inputReaderT1->SetFileName(T1FileName);
+ inputImageFilenames.push_back(T1FileName);
+ inputReaderT1->Update();
+ typename RescaleImageFilterType::Pointer rescalerT1 = RescaleImageFilterType::New();
+ // Set the upper limit to 4096 in the case of floating point data.
+ const inputPixelType outMin = 0;
+ const inputPixelType outMax
+ = ( vcl_numeric_limits<inputPixelType>::max() > 4096 ) ? 4096 : vcl_numeric_limits<inputPixelType>::max();
+ rescalerT1->SetOutputMinimum(outMin);
+ rescalerT1->SetOutputMaximum(outMax);
+ rescalerT1->SetInput( inputReaderT1->GetOutput() );
+ rescalerT1->Update();
+ images.push_back( rescalerT1->GetOutput() );
+ }
+
+ if( T2FileName == "" )
+ {
+ std::cerr << " T2 file is missing.\n ";
+ }
+ else
+ {
+ typename InputImageReaderType::Pointer inputReaderT2 = InputImageReaderType::New();
+ inputReaderT2->SetFileName(T2FileName);
+ inputImageFilenames.push_back(T2FileName);
+ inputReaderT2->Update();
+ typename RescaleImageFilterType::Pointer rescalerT2 = RescaleImageFilterType::New();
+ // Set the upper limit to 4096 in the case of floating point data.
+ const inputPixelType outMin = 0;
+ const inputPixelType outMax
+ = ( vcl_numeric_limits<inputPixelType>::max() > 4096 ) ? 4096 : vcl_numeric_limits<inputPixelType>::max();
+ rescalerT2->SetOutputMinimum(outMin);
+ rescalerT2->SetOutputMaximum(outMax);
+ rescalerT2->SetInput( inputReaderT2->GetOutput() );
+ rescalerT2->Update();
+ images.push_back( rescalerT2->GetOutput() );
+ }
+
+ if( PDFileName == "" )
+ {
+ std::cerr << " PD file is missing.\n ";
+ }
+ else
+ {
+ typename InputImageReaderType::Pointer inputReaderPD = InputImageReaderType::New();
+ inputReaderPD->SetFileName(PDFileName);
+ inputImageFilenames.push_back(PDFileName);
+ inputReaderPD->Update();
+ typename RescaleImageFilterType::Pointer rescalerPD = RescaleImageFilterType::New();
+ // Set the upper limit to 4096 in the case of floating point data.
+ const inputPixelType outMin = 0;
+ const inputPixelType outMax
+ = ( vcl_numeric_limits<inputPixelType>::max() > 4096 ) ? 4096 : vcl_numeric_limits<inputPixelType>::max();
+ rescalerPD->SetOutputMinimum(outMin);
+ rescalerPD->SetOutputMaximum(outMax);
+ rescalerPD->SetInput( inputReaderPD->GetOutput() );
+ rescalerPD->Update();
+ images.push_back( rescalerPD->GetOutput() );
+ }
+ if( T1FileName == "" && T2FileName == "" && PDFileName == "" )
+ {
+ std::cerr << " We need at lease one or more target image.\n";
+ status = -1;
+ return status;
+ }
+ // Compute the brain outline in order to cut images later after bias
+ // correction
+ const unsigned int closingSize = 7;
+ const float otsuPercentileThreshold = 0.01;
+ // typename InputImageType::Pointer HeadOutlineMaskImage =
+ // FindLargestForgroundFilledMask<InputImageType>( images[0],
+ // otsuPercentileThreshold, closingSize );
+ typedef itk::LargestForegroundFilledMaskImageFilter<InputImageType> LFFMaskFilterType;
+ typename LFFMaskFilterType::Pointer LFF = LFFMaskFilterType::New();
+ LFF->SetInput(images[0]);
+ LFF->SetOtsuPercentileThreshold(otsuPercentileThreshold);
+ LFF->SetClosingSize(closingSize);
+ LFF->Update();
+ typename InputImageType::Pointer HeadOutlineMaskImage = LFF->GetOutput();
+
+ std::cerr << " Following image will be used as a reference target image \n";
+ std::cerr << " Input Image 1:: " << images[0]->GetLargestPossibleRegion().GetSize() << std::endl;
+
+ /** - read atlas images **/
+ /* atlas should have list of atlas for interest and background */
+ /* and those atlas should including prior weight. */
+ std::vector<std::string> PriorsList(priorsList);
+ std::vector<double> PriorWeights(priorsWeightList);
+
+ typedef itk::Image<priorPixelType, 3> PriorImageType;
+ typedef itk::ImageFileReader<PriorImageType> PriorImageReaderType;
+ typedef typename PriorImageType::Pointer PriorImagePointer;
+
+ std::vector<PriorImagePointer> priors;
+
+ if( PriorsList.size() != PriorWeights.size() )
+ {
+ std::cerr << " Number of List of Prior Image should be same to the Prior Weights List\n";
+ status = -1;
+ return status;
+ }
+
+ std::vector<std::string> priorImageFileNames;
+ if( PriorsList.empty() )
+ {
+ std::cerr << " There is no prior image at all, we need at least one image\n";
+ status = -1;
+ return status;
+ }
+ else
+ {
+ for( unsigned int i = 0; i < PriorsList.size(); i++ )
+ {
+ typename PriorImageReaderType::Pointer priorReader = PriorImageReaderType::New();
+ std::string name = PriorsList.at(i);
+ priorReader->SetFileName( name );
+ priorImageFileNames.push_back(name);
+ priorReader->Update();
+ priors.push_back( priorReader->GetOutput() );
+ std::cout << " Prior List :: " << i << " :: " << name << "\n";
+ }
+ }
+
+ typedef EMSegmentationFilter<InputImageType, PriorImageType> SegFilterType;
+ typename SegFilterType::VectorType priorweights( PriorWeights.size() );
+ for( unsigned int i = 0; i < PriorWeights.size(); i++ )
+ {
+ priorweights[i] = PriorWeights.at(i);
+ }
+
+ std::cerr << "Start segmentation...\n";
+
+ typename SegFilterType::Pointer segfilter = SegFilterType::New();
+ segfilter->DebugOn();
+
+ // HACK: Need for loop around this
+ std::vector<InputImagePointer> templateImageList;
+ templateImageList.clear();
+
+ templateImageList.push_back( templateReader->GetOutput() );
+
+ segfilter->SetTemplateImages( templateImageList );
+ segfilter->SetInputImages( images );
+ segfilter->SetPriors( priors );
+ segfilter->SetMaximumIterations( maxIteration );
+ segfilter->SetPriorWeights(priorweights);
+ segfilter->SetMaxBiasDegree(degreeOfBiasFieldCorrection);
+ if( likelihoodTolerance > 0.0 ) // override the default only if provided.
+ {
+ segfilter->SetLikelihoodTolerance(likelihoodTolerance);
+ }
+ if( warp )
+ {
+ segfilter->DoWarpingOn();
+ }
+ else
+ {
+ segfilter->DoWarpingOff();
+ }
+ segfilter->Update();
+
+ // Write the labels
+ std::cerr << "Writing labels...\n";
+ {
+ typedef itk::ImageFileWriter<ByteImageType> OutputWriterType;
+ OutputWriterType::Pointer writer = OutputWriterType::New();
+
+ writer->SetInput( segfilter->GetOutput() );
+ writer->UseCompressionOn();
+ std::string fn = std::string(OutputFileNamePrefix + "_labels.nii.gz");
+ writer->SetFileName( fn.c_str() );
+ writer->Update();
+ }
+
+ // Write the secondary outputs
+ if( true )
+ {
+ std::cerr << "Writing filtered and bias corrected images...\n";
+ std::vector<InputImagePointer> imgset = segfilter->GetCorrected();
+ for( unsigned i = 0; i < imgset.size(); i++ )
+ {
+ typedef itk::CastImageFilter<InputImageType, ShortImageType> CasterType;
+ typename CasterType::Pointer caster = CasterType::New();
+ typename itk::MultiplyImageFilter<InputImageType,
+ InputImageType>::Pointer multiplyFilter
+ = itk::MultiplyImageFilter<InputImageType, InputImageType>::New();
+ multiplyFilter->SetInput1(imgset[i]);
+ multiplyFilter->SetInput2(HeadOutlineMaskImage);
+
+ caster->SetInput( multiplyFilter->GetOutput() );
+ caster->Update();
+ std::string fn = std::string(inputImageFilenames[i] + "_corrected.nii.gz");
+
+ typedef itk::ImageFileWriter<ShortImageType> ShortWriterType;
+ ShortWriterType::Pointer writer = ShortWriterType::New();
+
+ writer->SetInput( caster->GetOutput() );
+ writer->SetFileName( fn.c_str() );
+ writer->UseCompressionOn();
+ writer->Update();
+ }
+
+ // Short posteriors
+ std::cerr << "Writing posterior images...\n";
+ std::vector<ShortImagePointer> probset = segfilter->GetShortPosteriors();
+ for( unsigned int i = 0; i < ( probset.size() - 3 ); i++ )
+ {
+ typedef itk::ImageFileWriter<ShortImageType> ShortWriterType;
+ ShortWriterType::Pointer writer = ShortWriterType::New();
+ writer->SetInput( probset[i] );
+ writer->SetFileName(priorImageFileNames[i] + "_posterior.nii.gz");
+ writer->UseCompressionOn();
+ writer->Update();
+ }
+ }
+ return status;
+}
+
+int main(int argc, char *argv[])
+{
+ PARSE_ARGS;
+ int status = -1;
+ status
+ = simpleRunEMS<float, float>(t1Volume,
+ t2Volume,
+ pdVolume,
+ templateVolume,
+ maxIteration,
+ OutputFileNamePrefix,
+ priorsList,
+ priorsWeightList,
+ warp,
+ degreeOfBiasFieldCorrection,
+ likelihoodTolerance );
+ /*
+ if(inputPixelType=="double") status =
+ simpleRunEMS<double>( t1Volume, t2Volume, pdVolume, templateVolume, maxIteration ,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="short") status =
+ simpleRunEMS<short>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="ushort") status =
+ simpleRunEMS<unsigned short>(t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="int") status =
+ simpleRunEMS<int>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="uint") status =
+ simpleRunEMS<unsigned int>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="long") status =
+ simpleRunEMS<long>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="float") status =
+ simpleRunEMS<float>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else if( inputPixelType=="uchar") status =
+ simpleRunEMS<unsigned char>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ else( inputPixelType=="char") status =
+ simpleRunEMS<char>( t1Volume, t2Volume, pdVolume, templateVolume,
+ OutputFileNamePrefix, priorsList, priorsWeightList, warp , degreeOfBiasFieldCorrection );
+ */
+ return status;
+}
+
View
143 BRAINSABC/TestSuite/simpleEM.xml
@@ -0,0 +1,143 @@
+<?xml version="1.0" encoding="utf-8"?>
+<executable>
+ <category>UNC</category>
+ <title>simpleEM</title>
+ <description>Atlas-based tissue segmentation method </description>
+
+<!-- ## testing for simpleEM -->
+ <parameters>
+ <label>Input Images</label>
+ <image>
+ <name>t1Volume</name>
+ <longflag>t1Volume</longflag>
+ <description>T1 Volume</description>
+ <label>T1 Volume</label>
+ <default></default>
+ <channel>input</channel>
+ </image>
+ <image>
+ <name>t2Volume</name>
+ <longflag>t2Volume</longflag>
+ <description>T2 Volume</description>
+ <label>T2 Volume</label>
+ <default></default>
+ <channel>input</channel>
+ </image>
+ <image>
+ <name>pdVolume</name>
+ <longflag>pdVolume</longflag>
+ <description>PD Volume</description>
+ <label>PD Volume</label>
+ <default></default>
+ <channel>input</channel>
+ </image>
+ </parameters>
+ <parameters>
+ <label>Atlas Images</label>
+ <image>
+ <name>templateVolume</name>
+ <longflag>templateVolume</longflag>
+ <description>template Volume</description>
+ <label>template Volume</label>
+ <default></default>
+ <channel>input</channel>
+ </image>
+ <string-vector>
+ <name>priorsList</name>
+ <longflag>priorsList</longflag>
+ <description> number of input should same to the number of priorsList</description>
+ <default></default>
+ </string-vector>
+ <double-vector>
+ <name>priorsWeightList</name>
+ <longflag deprecatedalias="priorWeightsList">priorsWeightList</longflag>
+ <description> number of input should same to the number of priorsList</description>
+ <default></default>
+ </double-vector>
+ <double>
+ <name>likelihoodTolerance</name>
+ <longflag>likelihoodTolerance</longflag>
+ <description> a convergence parameter (0 indicates take the default from the EMSegmentationFilter constructor.)</description>
+ <default>0.0</default>
+ </double>
+ </parameters>
+ <parameters>
+ <label>Algorithm Parameters</label>
+ <boolean>
+ <name>warp</name>
+ <longflag>warp</longflag>
+ <label>warping </label>
+ <default>false</default>
+ <channel>input</channel>
+ </boolean>
+ <integer>
+ <name>degreeOfBiasFieldCorrection</name>
+ <longflag>degreeOfBiasFieldCorrection</longflag>
+ <label>maximum degree Of BiasField Correction</label>
+ <default>4</default>
+ <channel>input</channel>
+ </integer>
+ <integer>
+ <name>maxIteration</name>
+ <longflag>maxIteration</longflag>
+ <label>maximum iteration number</label>
+ <default>100</default>
+ <channel>input</channel>
+ </integer>
+ </parameters>
+ <parameters>
+ <label>IO data formats</label>
+ <string>
+ <longflag>OutputFileNamePrefix</longflag>
+ <label>Output File Prefix</label>
+ <default>BRAINSABC</default>
+ <channel>input</channel>
+ </string>
+ <string-enumeration>
+ <name>inputPixelType</name>
+ <longflag>inputPixelType</longflag>
+ <label>Input Type</label>
+ <description>Input Type</description>
+ <default>uchar</default>
+ <element>usigned char</element>
+ <element>short</element>
+ <element>unsigned short</element>
+ <element>int</element>
+ <element>unsigned int</element>
+ <element>long</element>
+ <element>float</element>
+ <element>double</element>
+ </string-enumeration>
+ <string-enumeration>
+ <name>outputPixelType</name>
+ <longflag>outputPixelType</longflag>
+ <label>Output Pixel Type</label>
+ <description>Input Type</description>
+ <default>uchar</default>
+ <element>unsigbed char</element>
+ <element>short</element>
+ <element>unsigned short</element>
+ <element>int</element>
+ <element>unsigned int</element>
+ <element>long</element>
+ <element>float</element>
+ <element>double</element>
+ </string-enumeration>
+ <string-enumeration>
+ <name>priorPixelType</name>
+ <longflag>priorPixelType</longflag>
+ <label>Prior Pixel Type</label>
+ <description>prior Type</description>
+ <default>float</default>
+ <element>unsigbed char</element>
+ <element>short</element>
+ <element>unsigned short</element>
+ <element>int</element>
+ <element>unsigned int</element>
+ <element>long</element>
+ <element>float</element>
+ <element>double</element>
+ </string-enumeration>
+ </parameters>
+ <!-- ## end of simpleEM -->
+</executable>
View
1  BRAINSABC/TestSuite/test_data/000517510_000517510_lmks.fcsv.md5
@@ -0,0 +1 @@
+91e369151012a32934b6218189a6d15d
View
1  BRAINSABC/TestSuite/test_data/EMS-Param.xml.md5
@@ -0,0 +1 @@
+f2d6f5a91deed770db64bc7204ec4be2
View
1  BRAINSABC/TestSuite/test_data/EMS.log.md5
@@ -0,0 +1 @@
+b6a47cf38e7e46f9f8496155546f0d60
View
1  BRAINSABC/TestSuite/test_data/EMS.xml.md5
@@ -0,0 +1 @@
+2a39f7a4ce55aecb1c491f54ad9ad493
View
1  BRAINSABC/TestSuite/test_data/FAST_logger.log.md5
@@ -0,0 +1 @@
+aa78504cd5f2401582b90346e42502c3
View
1  BRAINSABC/TestSuite/test_data/ISO_T1_REP0.nii.gz.md5
@@ -0,0 +1 @@
+67793308fc0d9fe80a6baf917a60cfb4
View
1  BRAINSABC/TestSuite/test_data/ISO_T1_REP0_to_template_EMS.bspline.md5
@@ -0,0 +1 @@
+927ddcdf4ab6048c041c4610d00789be
View
1  BRAINSABC/TestSuite/test_data/ISO_T1_REP1.nii.gz.md5
@@ -0,0 +1 @@
+9f074e4025a2054bf85a5d0c23d9ff2f
View
1  BRAINSABC/TestSuite/test_data/ISO_T2_REP0.nii.gz.md5
@@ -0,0 +1 @@
+d66e002dedb39db3062f364c7cf84cd4
View
1  BRAINSABC/TestSuite/test_data/T1_REP0.nii.gz.md5
@@ -0,0 +1 @@
+64498adb2f9d82743d4e4362acf4c08c
View
1  BRAINSABC/TestSuite/test_data/T1_REP1.nii.gz.md5
@@ -0,0 +1 @@
+2ae8d425cb4b7f522b915310648283d6
View
1  BRAINSABC/TestSuite/test_data/T2_REP0.nii.gz.md5
@@ -0,0 +1 @@
+ea059316f49ff3c0dbb4fef52b4aab26
View
1  BRAINSABC/TestSuite/test_data/binary_release.test.md5
@@ -0,0 +1 @@
+86415f5e364e7d2b87a0b4b9d7d6b9c9
View
1  BRAINSABC/TestSuite/test_data/medium_T1_REP0.nii.gz.md5
@@ -0,0 +1 @@
+d52bfda9686a0d83d88339dcb3fbe4da
View
1  BRAINSABC/TestSuite/test_data/medium_T1_REP1.nii.gz.md5
@@ -0,0 +1 @@
+f304199d85ca408a43db9bbaba0481cd
View
1  BRAINSABC/TestSuite/test_data/medium_T2_REP0.nii.gz.md5
@@ -0,0 +1 @@
+5610c3b62db01b09e252415c148a067e
View
1  BRAINSABC/TestSuite/test_data/qqq.logger.md5
@@ -0,0 +1 @@
+629c475ae1e6cad9b08b50a6c071ea63
View
1  BRAINSABC/TestSuite/test_data/refResults/labels.nii.gz.md5
@@ -0,0 +1 @@
+cc3f238174209038d79a076ac463b555
View
1  BRAINSABC/TestSuite/test_data/refResults/medium_T1_REP0_labels_BRAINSABC.nii.gz.md5
@@ -0,0 +1 @@
+2ddd57eba365997b63d69da79f7b15f1
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_corrected_BRAINSABC.nii.gz.md5
@@ -0,0 +1 @@
+f12d520b8806800584163b5bc8e7cb64
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_corrected_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+b106cbbd494a5d5bf6ce15967342cde0
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_labels_BRAINSABC.nii.gz.md5
@@ -0,0 +1 @@
+822f9b0f90bda6605933ace6502c8a4a
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_labels_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+37f384e7bb117dc5fef7527440da40d4
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_posterior0_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+4eeb4bf5e5bb00072f5a6d8b2686f03b
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_posterior1_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+7da5c82116eeceb5e7b290732c8b6e69
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_posterior2_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+79d40e6226122822d848683bea798ca1
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_registered_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+aa81874f49d7e74b229f9729959f800f
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_template_affine_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+e2c34ce16c5f801f6d3620a59c4fd140
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP0_template_warped_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+40838092d23ef6f67445aee701679048
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP1_corrected_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+65da92f6debbdb810658ac485ffc9553
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T1_REP1_registered_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+43a6cdbeff841b1f69db1c97ea90555a
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T2_REP0_corrected_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+03348131aa6626e2d6a795c433efa439
View
1  BRAINSABC/TestSuite/test_data/refResults/small_ISO_T2_REP0_registered_BRAINSABC.nrrd.md5
@@ -0,0 +1 @@
+51ba4d6ed6c93a7c6a80d0b431febf92
View
1  BRAINSABC/TestSuite/test_data/run_test.sh.md5
@@ -0,0 +1 @@
+22ea386f7bc25e4c623631ac3e9f8cd4
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0.nii.gz.md5
@@ -0,0 +1 @@
+f5376fed7ace5288f354a400f78b63c1
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_small_ISO_T1_REP1_EMS.affine.md5
@@ -0,0 +1 @@
+bcf7ac034b51a40c2a566074f7d7ac0d
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_small_ISO_T2_REP0_EMS.affine.md5
@@ -0,0 +1 @@
+e99fc06cd42c0f718b1ebac784e5a6c8
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_template_EMS.affine.md5
@@ -0,0 +1 @@
+aca942f1ae529e13b55ec0abe8c3809c
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP0_to_template_EMS.bspline.md5
@@ -0,0 +1 @@
+85f0d61e57edcf3e8632ad6a94c2059f
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T1_REP1.nii.gz.md5
@@ -0,0 +1 @@
+1efd6c51ce6846d04cd882778ae5982c
View
1  BRAINSABC/TestSuite/test_data/small_ISO_T2_REP0.nii.gz.md5
@@ -0,0 +1 @@
+29e6f3bda04965df586bdd37e16c3979
View
1  BRAINSABC/TestSuite/test_data/small_T1_REP0.nii.gz.md5
@@ -0,0 +1 @@
+5597e9dc203ffe689de3c0cc50cfee78
View
1  BRAINSABC/TestSuite/test_data/small_T1_REP0_to_template_EMS.bspline.md5
@@ -0,0 +1 @@
+1b5da6351f561587eb7e40edc6929952
View
1  BRAINSABC/TestSuite/test_data/small_T1_REP1.nii.gz.md5
@@ -0,0 +1 @@
+e94dcbc50400f911f328c218b4de7c83
View
1  BRAINSABC/TestSuite/test_data/small_T2_REP0.nii.gz.md5
@@ -0,0 +1 @@
+b7099ae59aedc0b7e5a9411d0e0368b1
View
1  BRAINSABC/TestSuite/test_data/test.sh.md5
@@ -0,0 +1 @@
+002e8928393b89364e40fc9b515f9b99
View
1  BRAINSABC/TestSuite/test_data/tmp_load.b2.md5
@@ -0,0 +1 @@
+99d28a03ddc239d339d148ba7a41e2ed
View
122 BRAINSABC/brainseg/AtlasCropImageSource.h
@@ -0,0 +1,122 @@
+//
+//
+// //////////////////////////////////////////////////////////////////////////////
+//
+// Generate cropped images using ROI obtained from probabilities
+//
+//
+//
+// //////////////////////////////////////////////////////////////////////////////
+
+// prastawa@cs.unc.edu 11/2003
+
+#ifndef __AtlasCropImageSource_h
+#define __AtlasCropImageSource_h
+
+#include "itkObject.h"
+
+#include <vector>
+
+/** \class AtlasCropImageSource
+ */
+template <class TInputImage, class TProbabilityImage>
+class AtlasCropImageSource : public itk::Object
+{
+public:
+
+ /** Standard class typedefs. */
+ typedef AtlasCropImageSource Self;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** The dimension of the image. */
+ itkStaticConstMacro(ImageDimension, unsigned int,
+ TInputImage::ImageDimension);
+
+ // Image types
+ typedef TInputImage InputImageType;
+ typedef typename TInputImage::Pointer InputImagePointer;
+ typedef typename TInputImage::IndexType InputImageIndexType;
+ typedef typename TInputImage::OffsetType InputImageOffsetType;
+ typedef typename TInputImage::PixelType InputImagePixelType;
+ typedef typename TInputImage::PointType InputImagePointType;
+ typedef typename TInputImage::RegionType InputImageRegionType;
+ typedef typename TInputImage::SizeType InputImageSizeType;
+ typedef typename TInputImage::SpacingType InputImageSpacingType;
+
+ typedef TProbabilityImage ProbabilityImageType;
+ typedef typename ProbabilityImageType::Pointer ProbabilityImagePointer;
+ typedef typename ProbabilityImageType::IndexType ProbabilityImageIndexType;
+ typedef typename ProbabilityImageType::OffsetType ProbabilityImageOffsetType;
+ typedef typename ProbabilityImageType::PixelType ProbabilityImagePixelType;
+ typedef typename ProbabilityImageType::RegionType ProbabilityImageRegionType;
+ typedef typename ProbabilityImageType::SizeType ProbabilityImageSizeType;
+ typedef typename ProbabilityImageType::SpacingType ProbabilityImageSpacingType;
+
+ typedef std::vector<ProbabilityImagePointer> ProbabilityImageList;
+
+ typedef struct
+ {
+ InputImageIndexType offset;
+ InputImageSizeType cropped_size;
+ InputImageSizeType original_size;
+ } CropInfoType;
+
+ // Set/get output image padding, in mm
+ itkGetMacro(Padding, double);
+ itkSetMacro(Padding, double);
+
+ bool CheckBounds();
+
+ void UseProbabilities(ProbabilityImageList probs);
+
+ // Create new images (either cropped or padded)
+ InputImagePointer Restore(InputImagePointer img);
+
+ InputImagePointer Crop(InputImagePointer img);
+
+ // Crop region information
+ void SetCropInfo(const CropInfoType & info)
+ {
+ m_CropInfo = info;
+ }
+
+ const CropInfoType & GetCropInfo()
+ {
+ return m_CropInfo;
+ }
+
+ // For debugging, generate slabs in last dim with top and bottom parts removed
+ itkSetMacro(SlabMode, bool);
+ itkGetConstMacro(SlabMode, bool);
+ itkBooleanMacro(SlabMode);
+protected:
+
+ AtlasCropImageSource();
+ ~AtlasCropImageSource()
+ {
+ }
+
+ double m_Padding;
+
+ InputImageIndexType m_LowerBound;
+ InputImageIndexType m_UpperBound;
+
+ InputImagePointType m_InputOrigin;
+ InputImagePointType m_CropOrigin;
+
+ InputImageSizeType m_OriginalSize;
+
+ CropInfoType m_CropInfo;
+
+ bool m_SlabMode;
+};
+
+#ifndef MU_MANUAL_INSTANTIATION
+#include "AtlasCropImageSource.txx"
+#endif
+
+#endif
View
291 BRAINSABC/brainseg/AtlasCropImageSource.txx
@@ -0,0 +1,291 @@
+#ifndef __AtlasCropImageSource_txx
+#define __AtlasCropImageSource_txx
+
+#include "itkImageRegionIteratorWithIndex.h"
+
+#include "AtlasCropImageSource.h"
+
+template <class TInputImage, class TProbabilityImage>
+AtlasCropImageSource<TInputImage, TProbabilityImage>
+::AtlasCropImageSource()
+{
+ m_Padding = 8.0;
+
+ m_LowerBound.Fill(0);
+ m_UpperBound.Fill(0);
+
+ m_OriginalSize.Fill(0);
+
+ m_SlabMode = false;
+}
+
+template <class TInputImage, class TProbabilityImage>
+bool
+AtlasCropImageSource<TInputImage, TProbabilityImage>
+::CheckBounds()
+{
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ if( m_LowerBound[i] > m_UpperBound[i] )
+ {
+ return false;
+ }
+ }
+
+ InputImageSizeType croppedSize;
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ croppedSize[i] = m_UpperBound[i] - m_LowerBound[i] + 1;
+ }
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ if( croppedSize[i] > m_OriginalSize[i] )
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+template <class TInputImage, class TProbabilityImage>
+void
+AtlasCropImageSource<TInputImage, TProbabilityImage>
+::UseProbabilities(ProbabilityImageList probs)
+{
+ if( probs.size() == 0 )
+ {
+ itkExceptionMacro(<< "Need at least one class probability image");
+ }
+
+ ProbabilityImageSizeType size
+ = probs[0]->GetLargestPossibleRegion().size();
+
+ ProbabilityImageSpacingType spacing
+ = probs[0]->GetSpacing();
+ // Make sure all class probabilities have the same space
+ for( unsigned int i = 1; i < probs.size(); i++ )
+ {
+ ProbabilityImageSizeType othersize
+ = probs[i]->GetLargestPossibleRegion().size();
+ if( size != othersize )
+ {
+ itkExceptionMacro(<< "Probability list size mismatch");
+ }
+ }
+
+ // Transform padding to voxel counts
+ InputImageOffsetType padding;
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ padding[i] = (unsigned int)vcl_floor(m_Padding / spacing[i] + 0.5);
+ }
+ // Make sure padding is sensible
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ if( (long)size[i] <= padding[i] )
+ {
+ itkExceptionMacro(
+ << "Bounding box padding larger than or equal to image size");
+ }
+ }
+
+ // Save size info
+ m_OriginalSize = size;
+ // Initial bounds: whole image
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ m_LowerBound[i] = size[i] - 1;
+ m_UpperBound[i] = 0;
+ }
+
+ // Go through image and update bounds
+ typedef itk::ImageRegionIteratorWithIndex<ProbabilityImageType>
+ IteratorType;
+
+ IteratorType it( probs[0], probs[0]->GetLargestPossibleRegion() );
+
+ it.GoToBegin();
+ while( !it.IsAtEnd() )
+ {
+ ProbabilityImageIndexType ind = it.GetIndex();
+
+ double sumProb = 0;
+ for( unsigned int i = 0; i < probs.size(); i++ )
+ {
+ sumProb += probs[i]->GetPixel(ind);
+ }
+
+ if( sumProb > 0 )
+ {
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ if( ind[i] < m_LowerBound[i] )
+ {
+ m_LowerBound[i] = ind[i];
+ }
+ if( ind[i] > m_UpperBound[i] )
+ {
+ m_UpperBound[i] = ind[i];
+ }
+ }
+ }
+
+ ++it;
+ }
+
+ // Only generate a slab?
+ if( m_SlabMode )
+ {
+ long len
+ = m_UpperBound[ImageDimension - 1] - m_LowerBound[ImageDimension - 1] + 1;
+ long offt = (long)( 0.2 * len );
+ m_LowerBound[ImageDimension - 1] += offt;
+ m_UpperBound[ImageDimension - 1] -= offt;
+ }
+ // Enlarge/pad bounding box
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ if( m_LowerBound[i] > padding[i] )
+ {
+ m_LowerBound[i] -= padding[i];
+ }
+ else
+ {
+ m_LowerBound[i] = 0;
+ }
+
+ if( m_UpperBound[i] < ( (long)size[i] - 1 - padding[i] ) )
+ {
+ m_UpperBound[i] += padding[i];
+ }
+ else
+ {
+ m_UpperBound[i] = size[i] - 1;
+ }
+ }
+
+ // Store origin information
+ m_InputOrigin = probs[0]->GetOrigin();
+ probs[0]->TransformIndexToPhysicalPoint(m_LowerBound, m_CropOrigin);
+
+ // m_CropInfo.offset =
+ // m_CropInfo.size =
+}
+
+template <class TInputImage, class TProbabilityImage>
+typename AtlasCropImageSource<TInputImage, TProbabilityImage>
+::InputImagePointer
+AtlasCropImageSource<TInputImage, TProbabilityImage>
+::Restore(InputImagePointer img)
+{
+ if( !this->CheckBounds() )
+ {
+ itkExceptionMacro(<< "Invalid bounds");
+ }
+
+ InputImageSizeType size = img->GetLargestPossibleRegion().size();
+
+ InputImageSizeType croppedSize;
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ croppedSize[i] = m_UpperBound[i] - m_LowerBound[i] + 1;
+ }
+
+ if( size != croppedSize )
+ {
+ itkExceptionMacro(<< "Input size does not match size of cropped image");
+ }
+
+ InputImageRegionType region;
+ region.SetSize(m_OriginalSize);
+
+ InputImagePointer output = InputImageType::New();
+ output->CopyInformation(img);
+ // Adjust origin
+ output->SetOrigin(m_InputOrigin);
+ output->SetRegions(region);
+ output->Allocate();
+ output->FillBuffer(0);
+
+ InputImageOffsetType offt;
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ offt[i] = m_LowerBound[i];
+ }
+
+ typedef itk::ImageRegionIteratorWithIndex<InputImageType> IteratorType;
+
+ IteratorType inIt( img, img->GetLargestPossibleRegion() );
+
+ inIt.GoToBegin();
+ while( !inIt.IsAtEnd() )
+ {
+ InputImageIndexType ind = inIt.GetIndex();
+
+ output->SetPixel( ind + offt, inIt.Get() );
+
+ ++inIt;
+ }
+
+ return output;
+}
+
+template <class TInputImage, class TProbabilityImage>
+typename AtlasCropImageSource<TInputImage, TProbabilityImage>
+::InputImagePointer
+AtlasCropImageSource<TInputImage, TProbabilityImage>
+::Crop(InputImagePointer img)
+{
+ if( !this->CheckBounds() )
+ {
+ itkExceptionMacro(<< "Invalid bounds");
+ }
+
+ InputImageSizeType size = img->GetLargestPossibleRegion().size();
+
+ if( size != m_OriginalSize )
+ {
+ itkExceptionMacro(<< "Input size does not match size of probability images");
+ }
+
+ InputImageSizeType croppedSize;
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ croppedSize[i] = m_UpperBound[i] - m_LowerBound[i] + 1;
+ }
+
+ InputImageRegionType cropRegion;
+ cropRegion.SetSize(croppedSize);
+
+ InputImagePointer output = InputImageType::New();
+ output->CopyInformation(img);
+ // Adjust origin
+ output->SetOrigin(m_CropOrigin);
+ output->SetRegions(cropRegion);
+ output->Allocate();
+
+ InputImageOffsetType offt;
+ for( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ offt[i] = m_LowerBound[i];
+ }
+
+ typedef itk::ImageRegionIteratorWithIndex<InputImageType> IteratorType;
+
+ IteratorType outIt(output, cropRegion);
+
+ outIt.GoToBegin();
+ while( !outIt.IsAtEnd() )
+ {
+ InputImageIndexType ind = outIt.GetIndex();
+
+ outIt.Set( img->GetPixel(ind + offt) );
+
+ ++outIt;
+ }
+
+ return output;
+}
+
+#endif
View
384 BRAINSABC/brainseg/AtlasDefinition.cxx
@@ -0,0 +1,384 @@
+#include "AtlasDefinition.h"
+#include <expat.h>
+#include <iostream>
+#include <fstream>
+#include <stdlib.h>
+
+namespace // anon
+{
+extern "C"
+{
+void
+start(void *data, const char *el, const char * *)
+{
+ // std::cerr << "Start, El = ("
+ // << el
+ // << ")";
+ // std::cerr << " Attributes = (";
+ // for(unsigned i = 0; attr[i] != 0; i+= 2)
+ // {
+ // if(i > 0) std::cerr << ",";
+ // std::cerr << attr[i] << "=" << attr[i+1];
+ // }
+ // std::cerr << ")" << std::endl;
+ AtlasDefinition *_this = reinterpret_cast<AtlasDefinition *>( data );
+
+ _this->XMLStart(el);
+}
+
+void
+end(void *data, const char *el)
+{
+ // std::cerr << "End, El = ("
+ // << el
+ // << ")" << std::endl;
+ AtlasDefinition *_this = reinterpret_cast<AtlasDefinition *>( data );
+
+ _this->XMLEnd(el);
+}
+
+void
+charhandler(void *data, const char *txt, int txtlen)
+{
+ // std::cerr << "Char data = (";
+ // for(unsigned i = 0; i < txtlen; i++)
+ // {
+ // std::cerr << txt[i];
+ // }
+ // std::cerr << ")" << std::endl;
+ char *buf = new char[txtlen + 1];
+ int i;
+ for( i = 0; i < txtlen; i++ )
+ {
+ buf[i] = txt[i];
+ }
+ buf[i] = '\0';
+ AtlasDefinition *_this = reinterpret_cast<AtlasDefinition *>( data );
+ _this->XMLChar(buf);
+ delete[] buf;
+}
+
+}
+}
+#if 0
+const char *AtlasDefinition::tissueTypes[] =
+ {
+ "WM",
+ "GM",
+ "BGM",
+ "CSF",
+ "VB",
+ "NOTCSF",
+ "NOTGM",
+ "NOTVB",
+ "NOTWM",
+ "AIR",
+ 0
+ };
+#endif
+
+AtlasDefinition::AtlasDefinition() :
+ m_LastWeight(0.0),
+ m_LastLower(0.0),
+ m_LastUpper(0.0),
+ m_LastGaussianClusterCount(0),
+ m_LastLabelCode(0),
+ m_LastUseForBias(0),
+ m_LastIsForegroundPrior(0)
+{
+ this->m_TissueTypes.resize(0);
+}
+
+const AtlasDefinition::TissueTypeVector &
+AtlasDefinition::TissueTypes()
+{
+ if( this->m_TissueTypes.size() == 0 )
+ {
+ for( PriorMapType::const_iterator it = this->m_PriorMap.begin();
+ it != this->m_PriorMap.end();
+ ++it )
+ {
+ this->m_TissueTypes.push_back(it->first);
+ }
+ }
+ return this->m_TissueTypes;
+}
+
+void
+AtlasDefinition::XMLStart(const char *el)
+{
+ std::string El(el);
+
+ this->m_XMLElementStack.push_back(El);
+ if( El == "Prior" )
+ {
+ this->m_LastPriorBounds.clear();
+ }
+}
+
+void
+AtlasDefinition::XMLEnd(const char *el)
+{
+ std::string El(el);
+ const char *start = this->m_LastXMLString.c_str();
+ char * last;
+
+ // pop the current element name off the stack.
+ this->m_XMLElementStack.pop_back();
+
+ if( El == "Atlas" )
+ {
+ return; // final element
+ }
+ const std::string ContainingElement = this->m_XMLElementStack.back();
+ if( El == "filename" )
+ {
+ this->m_LastFilename = this->m_LastXMLString;
+ }
+ else if( El == "type" )
+ {
+ if( ContainingElement == "Prior" )
+ {
+ this->m_LastPriorType = this->m_LastXMLString;
+ }
+ else
+ {
+ this->m_LastType = this->m_LastXMLString;
+ }
+ }
+ else if( El == "AtlasImage" )
+ {
+ this->m_TemplateVolumes[this->m_LastType] = this->m_LastFilename;
+ }
+ else if( El == "Weight" )
+ {
+ double w = strtod(start, &last);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad Weight given " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastWeight = w;
+ }
+ else if( El == "lower" )
+ {
+ double l = strtod(start, &last);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad lower bound " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastLower = l;
+ }
+ else if( El == "upper" )
+ {
+ double u = strtod(start, &last);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad upper bound " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastUpper = u;
+ }
+ else if( El == "GaussianClusterCount" )
+ {
+ int GCC = strtol(start, &last, 10);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad GaussianClusterCount given " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastGaussianClusterCount = GCC;
+ }
+ else if( El == "LabelCode" )
+ {
+ int GCC = strtol(start, &last, 10);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad LabelCode given " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastLabelCode = GCC;
+ }
+ else if( El == "UseForBias" )
+ {
+ int UFB = strtol(start, &last, 10);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad UseForBias given " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastUseForBias = UFB;
+ }
+ else if( El == "IsForegroundPrior" )
+ {
+ int UFB = strtol(start, &last, 10);
+ if( start == (const char *)last )
+ {
+ std::cerr << "Bad IsForegroundPrior given " << this->m_LastXMLString
+ << std::endl;
+ throw;
+ }
+ this->m_LastIsForegroundPrior = UFB;
+ }
+ else if( El == "BrainMask" )
+ {
+ this->m_TemplateBrainMask = this->m_LastFilename;
+ }
+ else if( El == "HeadRegion" )
+ {
+ this->m_TemplateHeadRegion = this->m_LastFilename;
+ }
+ else if( El == "Prior" )
+ {
+ Prior curPrior;
+ curPrior.SetFilename(this->m_LastFilename);
+ curPrior.SetWeight(this->m_LastWeight);
+ curPrior.SetGaussianClusterCount(this->m_LastGaussianClusterCount);
+ curPrior.SetLabelCode(this->m_LastLabelCode);
+ curPrior.SetUseForBias(this->m_LastUseForBias);
+ curPrior.SetIsForegroundPrior(this->m_LastIsForegroundPrior);
+ curPrior.SetBoundsList(this->m_LastPriorBounds);
+ this->m_PriorMap[this->m_LastPriorType] = curPrior;
+ }
+ else if( El == "bounds" )
+ {
+ BoundsType & curBounds = this->m_LastPriorBounds[this->m_LastType];
+ curBounds.SetLower(this->m_LastLower);
+ curBounds.SetUpper(this->m_LastUpper);
+ }
+ else
+ {
+ std::cerr << "Unhandled XML Element type " << El << std::endl;
+ throw;
+ }
+}
+
+void
+AtlasDefinition::XMLChar(const char *buf)
+{
+ this->m_LastXMLString = buf;
+}
+
+void
+AtlasDefinition::InitFromXML(const std::string & XMLFilename)
+{
+ std::ifstream xmlFile(XMLFilename.c_str(),
+ std::ifstream::in);
+
+ if( !xmlFile.good() )
+ {
+ throw;
+ }
+ xmlFile.seekg(static_cast<std::streampos>( 0 ),
+ std::ios_base::end);
+ std::streampos fSize( xmlFile.tellg() );
+ xmlFile.seekg(static_cast<std::streampos>( 0 ),
+ std::ios_base::beg);
+ XML_Parser parser = XML_ParserCreate(0);
+ XML_SetUserData( parser, static_cast<void *>( this ) );
+ XML_SetElementHandler(parser, start, end);
+ XML_SetCharacterDataHandler(parser, charhandler);
+ char *filebuf = new char[fSize];
+ if( filebuf == NULL )
+ {
+ throw;
+ }
+ xmlFile.read(filebuf, fSize);
+ if( !xmlFile.good() )
+ {
+ delete[] filebuf;
+ throw;
+ }
+ xmlFile.close();
+ if( XML_Parse(parser, filebuf, fSize, 1) == 0 )
+ {
+ delete[] filebuf;
+ std::cerr << "XML File parsing error" << std::endl;
+ throw;
+ }
+ delete[] filebuf;
+}
+
+void
+AtlasDefinition::DebugPrint()
+{
+ std::cout << "<Atlas>" << std::endl;
+ for( TemplateMap::iterator it
+ = m_TemplateVolumes.begin();
+ it != m_TemplateVolumes.end();
+ ++it )
+ {
+ std::cout << " <AtlasImage>" << std::endl
+ << " <type>" << it->first << "</type>" << std::endl
+ << " <filename>" << it->second << "</filename>" << std::endl
+ << " </AtlasImage>" << std::endl;
+ }
+
+ std::cout << " <BrainMask>" << std::endl
+ << " <filename>"
+ << m_TemplateBrainMask
+ << "</filename>" << std::endl
+ << " </BrainMask>" << std::endl;
+
+ std::cout << " <HeadRegion>" << std::endl
+ << " <filename>"
+ << m_TemplateHeadRegion
+ << "</filename>" << std::endl
+ << " </HeadRegion>" << std::endl;
+ for( PriorMapType::iterator it
+ = m_PriorMap.begin();
+ it != m_PriorMap.end();
+ ++it )
+ {
+ std::cout << " <Prior>" << std::endl
+ << " <type>"
+ << it->first
+ << "</type>" << std::endl
+ << " <filename>"
+ << it->second.GetFilename()
+ << "</filename>" << std::endl
+ << " <Weight>"
+ << it->second.GetWeight()
+ << "</Weight>" << std::endl
+ << " <GaussianClusterCount>"
+ << it->second.GetGaussianClusterCount()
+ << "</GaussianClusterCount>" << std::endl
+ << " <LabelCode>"
+ << it->second.GetLabelCode()
+ << "</LabelCode>" << std::endl
+ << " <UseForBias>"
+ << it->second.GetUseForBias()
+ << "</UseForBias>" << std::endl
+ << " <IsForegroundPrior>"
+ << it->second.GetIsForegroundPrior()
+ << "</IsForegroundPrior>" << std::endl;
+ const BoundsMapType & curBounds = it->second.GetBoundsList();
+ for( BoundsMapType::const_iterator it2
+ = curBounds.begin();
+ it2 != curBounds.end();
+ ++it2 )
+ {
+ std::cout << " <bounds>" << std::endl
+ << " <type>"
+ << it2->first
+ << "</type>" << std::endl
+ << " <lower>"
+ << it2->second.GetLower()
+ << "</lower>" << std::endl
+ << " <upper>"
+ << it2->second.GetUpper()
+ << "<upper>" << std::endl
+ << " </bounds>" << std::endl;
+ }
+ std::cout << " </Prior>" << std::endl;
+ }
+ std::cout << "</Atlas>" << std::endl;
+}
+
View
361 BRAINSABC/brainseg/AtlasDefinition.h
@@ -0,0 +1,361 @@
+#ifndef __AtlasDefinition_h
+#define __AtlasDefinition_h
+
+#include <string>
+#include <vector>
+#include <list>
+#include <map>
+#include <iostream>
+
+/** \class AtlasDefiniton */
+class AtlasDefinition
+{
+public:
+ // static const char * tissueTypes[];
+ typedef std::vector<std::string> TissueTypeVector;
+ typedef std::map<std::string, std::string> TemplateMap;
+ typedef std::vector<float> FloatVector;
+ typedef std::vector<int> IntVector;
+ /** \class BoundsType */
+ class BoundsType
+ {
+ public:
+ BoundsType() : m_Low(0), m_High(0)
+ {
+ }
+
+ void SetLower(const double v)
+ {
+ m_Low = v;
+ }
+
+ double GetLower() const
+ {
+ return m_Low;
+ }
+
+ void SetUpper(const double v)
+ {
+ m_High = v;
+ }
+
+ double GetUpper() const
+ {
+ return m_High;
+ }
+
+ void Print(void) const
+ {
+ std::cout << "RANGE: [" << m_Low << "," << m_High << "]" << std::endl;
+ }
+
+ private:
+ double m_Low;
+ double m_High;
+ };
+ typedef std::map<std::string, BoundsType> BoundsMapType;
+
+ AtlasDefinition();
+ void InitFromXML(const std::string & XMLFilename);
+
+ void DebugPrint();
+
+ const TemplateMap & GetTemplateVolumes() const
+ {
+ return m_TemplateVolumes;
+ }
+
+ const std::string & GetTemplateBrainMask() const
+ {
+ return m_TemplateBrainMask;
+ }
+
+ const std::string & GetTemplateHeadRegion() const
+ {
+ return m_TemplateHeadRegion;
+ }
+
+ std::string GetPriorFilename(const std::string & tissueType) const
+ {
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING TISSUE TYPE IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetFilename();
+ }
+
+ double GetWeight(const std::string & tissueType) const
+ {
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING TISSUE TYPE IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetWeight();
+ }
+
+ int GetGaussianClusterCount(const std::string & tissueType) const
+ {
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING TISSUE TYPE IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetGaussianClusterCount();
+ }
+
+ int GetLabelCode(const std::string & tissueType) const
+ {
+ //HACK: All the get functions need review to remove duplicate code.
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING LABEL CODE IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetLabelCode();
+ }
+
+ int GetUseForBias(const std::string & tissueType) const
+ {
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING TISSUE TYPE IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetUseForBias();
+ }
+
+ bool GetIsForegroundPrior(const std::string & tissueType) const
+ {
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING IsForegroungPrior IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetIsForegroundPrior();
+ }
+
+ const BoundsType & GetBounds(const std::string & tissueType,
+ const std::string & Modality) const
+ {
+ if( m_PriorMap.find(tissueType) == m_PriorMap.end() )
+ {
+ std::cout << "MISSING TISSUE TYPE IN ATLAS: " << tissueType << std::endl;
+ throw;
+ }
+ PriorMapType::const_iterator mit = m_PriorMap.find(tissueType);
+ if( mit == m_PriorMap.end() )
+ {
+ // HACK: Should throw and exception here with line number and file
+ std::cout << "ERROR: Invalid tissueType requested in GetPriorFilename" << std::endl;
+ throw;
+ }
+ return mit->second.GetBounds(Modality);
+ }
+
+ double GetLow(const std::string & tissueType,
+ const std::string & Modality) const
+ {
+ return this->GetBounds(tissueType, Modality).GetLower();
+ }
+
+ double GetHigh(const std::string & tissueType,
+ const std::string & Modality) const
+ {
+ return this->GetBounds(tissueType, Modality).GetUpper();
+ }
+
+ const TissueTypeVector & TissueTypes();
+
+ void XMLStart(const char *el);
+
+ void XMLEnd(const char *el);
+
+ void XMLChar(const char *buf);
+
+private:
+ TemplateMap m_TemplateVolumes;
+ std::string m_TemplateBrainMask;
+ std::string m_TemplateHeadRegion;
+ TissueTypeVector m_TissueTypes;
+ class Prior
+ {
+ public:
+ Prior() :
+ m_Weight(0.0),
+ m_GaussianClusterCount(0),
+ m_LabelCode(0),
+ m_UseForBias(false),
+ m_IsForegroundPrior(false)
+ {
+ }
+
+ const std::string & GetFilename() const
+ {
+ return m_Filename;
+ }
+
+ void SetFilename(const std::string & s)
+ {
+ m_Filename = s;
+ }
+
+ double GetWeight() const
+ {
+ return m_Weight;
+ }
+
+ void SetWeight(double x)
+ {
+ m_Weight = x;
+ }
+
+ int GetGaussianClusterCount() const
+ {
+ return m_GaussianClusterCount;
+ }
+
+ void SetGaussianClusterCount(int i)
+ {
+ m_GaussianClusterCount = i;
+ }
+
+ int GetLabelCode() const
+ {
+ return m_LabelCode;
+ }
+
+ void SetLabelCode(const int i)
+ {
+ m_LabelCode = i;
+ }
+
+ int GetUseForBias() const
+ {
+ return m_UseForBias;
+ }
+
+ void SetUseForBias(const bool i)
+ {
+ m_UseForBias = i;
+ }
+
+ int GetIsForegroundPrior() const
+ {
+ return m_IsForegroundPrior;
+ }
+
+ void SetIsForegroundPrior(const bool i)
+ {
+ m_IsForegroundPrior = i;
+ }
+
+ const BoundsType & GetBounds(const std::string & Modality) const
+ {
+ BoundsMapType::const_iterator bit= m_BoundsMap.find(Modality);
+ if(bit == m_BoundsMap.end() )
+ {
+ std::cout << "MISSING MODALIITY TYPE IN ATLAS: " << Modality << std::endl;
+ throw;
+ }
+
+ return bit->second;
+ }
+
+ void SetBounds(const std::string & Modality, const BoundsType & b)
+ {
+ this->m_BoundsMap[Modality] = b;
+ }
+
+ void SetBounds(const std::string & Modality, double lower, double upper)
+ {
+ BoundsType b;
+
+ b.SetLower(lower);
+ b.SetUpper(upper);
+ this->SetBounds(Modality, b);
+ }
+
+ void SetBoundsList(BoundsMapType & boundsMap)
+ {
+ this->m_BoundsMap = boundsMap;
+ }
+
+ const BoundsMapType & GetBoundsList() const
+ {
+ return this->m_BoundsMap;
+ }
+
+ private:
+ std::string m_Filename;
+ double m_Weight;
+ int m_GaussianClusterCount;
+ int m_LabelCode;
+ bool m_UseForBias;
+ bool m_IsForegroundPrior;
+ BoundsMapType m_BoundsMap;
+ };
+ typedef std::map<std::string, Prior> PriorMapType;
+ PriorMapType m_PriorMap;
+ //
+ // XML Parsing variables
+ std::string m_LastXMLString;
+ std::string m_LastPriorType;
+ std::string m_LastFilename;
+ std::string m_LastType;
+ std::list<std::string> m_XMLElementStack;
+ double m_LastWeight,
+ m_LastLower,
+ m_LastUpper;
+
+ int m_LastGaussianClusterCount;
+ int m_LastLabelCode;
+ bool m_LastUseForBias;
+ bool m_LastIsForegroundPrior;
+ BoundsMapType m_LastPriorBounds;
+};
+
+#endif // AtlasDefinition_h
View
195 BRAINSABC/brainseg/AtlasRegistrationMethod.h
@@ -0,0 +1,195 @@
+//
+//
+// //////////////////////////////////////////////////////////////////////////////
+//
+// Registration of a dataset to an atlas using affine transformation and
+// MI image match metric
+//
+// Only for 3D!
+//
+// Given a list of filenames for atlas template and probabilities along with
+// the dataset, this class generate images that are in the space of the first
+// image (all data and probability images).
+//
+//
+//
+// //////////////////////////////////////////////////////////////////////////////
+
+// prastawa@cs.unc.edu 10/2003
+
+#ifndef __AtlasRegistrationMethod_h
+#define __AtlasRegistrationMethod_h
+
+#include "itkAffineTransform.h"
+#include "itkArray.h"
+#include "itkImage.h"
+#include "itkObject.h"
+
+#include <vector>
+
+#include "itkAffineTransform.h"
+#include "BRAINSFitHelper.h"
+
+#include <string>
+
+/** \class AtlasRegistrationMethod
+ */
+template <class TOutputPixel, class TProbabilityPixel>
+class AtlasRegistrationMethod : public itk::Object
+{
+public:
+
+ /** Standard class typedefs. */
+ typedef AtlasRegistrationMethod Self;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ // Image types
+ typedef itk::Image<TOutputPixel, 3> OutputImageType;
+ typedef typename OutputImageType::Pointer OutputImagePointer;
+ typedef typename OutputImageType::IndexType OutputImageIndexType;
+ typedef typename OutputImageType::OffsetType OutputImageOffsetType;
+ typedef typename OutputImageType::PixelType OutputImagePixelType;
+ typedef typename OutputImageType::SizeType OutputImageSizeType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+ typedef itk::Image<TProbabilityPixel, 3> ProbabilityImageType;
+ typedef typename ProbabilityImageType::Pointer ProbabilityImagePointer;
+ typedef typename ProbabilityImageType::IndexType ProbabilityImageIndexType;
+ typedef typename ProbabilityImageType::OffsetType ProbabilityImageOffsetType;
+ typedef typename ProbabilityImageType::PixelType ProbabilityImagePixelType;
+ typedef typename ProbabilityImageType::SizeType ProbabilityImageSizeType;
+ typedef typename ProbabilityImageType::RegionType ProbabilityImageRegionType;
+
+ typedef itk::Image<float, 3> InternalImageType;
+ typedef typename InternalImageType::Pointer InternalImagePointer;
+ typedef typename InternalImageType::IndexType InternalImageIndexType;
+ typedef typename InternalImageType::OffsetType InternalImageOffsetType;
+ typedef typename InternalImageType::PixelType InternalImagePixelType;
+ typedef typename InternalImageType::RegionType InternalImageRegionType;
+ typedef typename InternalImageType::SizeType InternalImageSizeType;
+
+ typedef itk::Image<unsigned char, 3> ByteImageType;
+ typedef typename ByteImageType::Pointer ByteImagePointer;
+ typedef typename ByteImageType::IndexType ByteImageIndexType;
+ typedef typename ByteImageType::OffsetType ByteImageOffsetType;
+ typedef typename ByteImageType::PixelType ByteImagePixelType;
+ typedef typename ByteImageType::RegionType ByteImageRegionType;
+ typedef typename ByteImageType::SizeType ByteImageSizeType;
+
+ typedef std::vector<ProbabilityImagePointer> ProbabilityImageList;
+ typedef std::vector<OutputImagePointer> OutputImageList;
+
+ typedef itk::Array<unsigned char> FlagArrayType;
+
+ void SetSuffix(std::string suffix);
+
+ itkGetConstMacro(OutputDebugDir, std::string);
+ itkSetMacro(OutputDebugDir, std::string);
+
+ void SetAtlasOriginalImageList(std::vector<InternalImagePointer> & NewAtlasList);
+
+ void SetIntraSubjectOriginalImageList(std::vector<InternalImagePointer> & NewImageList);
+
+ // itkSetMacro( IntraSubjectTransformFileNames, std::vector<std::string> );
+ itkSetMacro( AtlasToSubjectTransformFileName, std::string );
+
+ void SetIntraSubjectTransformFileNames(std::vector<std::string> userlist)
+ {
+ m_IntraSubjectTransformFileNames = userlist;
+ }
+
+ void RegisterImages();
+
+ GenericTransformType::Pointer GetAtlasToSubjectTransform()
+ {
+ return m_AtlasToSubjectTransform;
+ }
+
+ std::vector<GenericTransformType::Pointer> GetIntraSubjectTransforms()
+ {
+ return m_IntraSubjectTransforms;
+ }
+
+ // Set/Get the Debugging level for filter verboseness
+ itkSetMacro(DebugLevel, unsigned int);
+ itkGetMacro(DebugLevel, unsigned int);
+
+ itkGetMacro(UseNonLinearInterpolation, bool);
+ itkSetMacro(UseNonLinearInterpolation, bool);
+
+ void Update();
+
+ void SetAtlasLinearTransformChoice(const std::string & c)
+ {
+ m_AtlasLinearTransformChoice = c;
+ }
+
+ void SetImageLinearTransformChoice(const std::string & c)
+ {
+ m_ImageLinearTransformChoice = c;
+ }
+
+ void SetWarpGrid(const unsigned int gx, const unsigned int gy, const unsigned int gz)
+ {
+ m_WarpGrid.resize(3); m_WarpGrid[0] = gx; m_WarpGrid[1] = gy; m_WarpGrid[2] = gz;
+ }
+
+ // Get and set the input volume image types. i.e. T1 or T2 or PD
+ void SetInputVolumeTypes(const std::vector<std::string> & newInputVolumeTypes)
+ {
+ this->m_InputVolumeTypes = newInputVolumeTypes;
+ }
+
+ std::vector<std::string> GetInputVolumeTypes(void) const
+ {
+ return this->m_InputVolumeTypes;
+ }
+
+protected:
+
+ AtlasRegistrationMethod();
+ ~AtlasRegistrationMethod();
+
+ OutputImagePointer CopyOutputImage(InternalImagePointer img);
+
+ ProbabilityImagePointer CopyProbabilityImage(InternalImagePointer img);
+
+private:
+
+ std::string m_Suffix;
+ std::string m_OutputDebugDir;
+
+ // ByteImagePointer m_AtlasOriginalMask;
+ std::vector<InternalImagePointer> m_AtlasOriginalImageList;
+ std::vector<InternalImagePointer> m_IntraSubjectOriginalImageList;
+ ByteImagePointer m_InputImageTissueRegion;
+ ImageMaskPointer m_InputSpatialObjectTissueRegion;
+
+ std::vector<unsigned int> m_WarpGrid;
+ std::vector<std::string> m_IntraSubjectTransformFileNames;
+ std::string m_AtlasToSubjectTransformFileName;
+
+ std::vector<std::string> m_InputVolumeTypes;
+
+ GenericTransformType::Pointer m_AtlasToSubjectTransform;
+ std::vector<GenericTransformType::Pointer> m_IntraSubjectTransforms;
+
+ bool m_UseNonLinearInterpolation;
+ bool m_DoneRegistration;
+ bool m_Modified;
+
+ std::string m_AtlasLinearTransformChoice;
+ std::string m_ImageLinearTransformChoice;
+
+ unsigned int m_DebugLevel;
+};
+
+#ifndef MU_MANUAL_INSTANTIATION
+#include "AtlasRegistrationMethod.txx"
+#endif
+
+#endif
View
772 BRAINSABC/brainseg/AtlasRegistrationMethod.txx
@@ -0,0 +1,772 @@
+#ifndef __AtlasRegistrationMethod_txx
+#define __AtlasRegistrationMethod_txx
+
+#include "itkAffineTransform.h"
+#include "itkBSplineInterpolateImageFunction.h"
+#include "itkImageFileReader.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "itkNearestNeighborInterpolateImageFunction.h"
+#include "itkResampleImageFilter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+// MI registration module
+#include "AtlasRegistrationMethod.h"
+#include "RegistrationParameters.h"
+
+#include "vnl/vnl_math.h"
+
+#include "Log.h"
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+
+#include "itkBRAINSROIAutoImageFilter.h"
+
+GenericTransformType::Pointer MakeRigidIdentity(void)
+{
+ // Also append identity matrix for each image
+ VersorRigid3DTransformType::Pointer rigidIdentity = VersorRigid3DTransformType::New();
+
+ rigidIdentity->SetIdentity();
+ GenericTransformType::Pointer genericTransform = NULL;
+ genericTransform = rigidIdentity.GetPointer();
+ return genericTransform;
+}
+
+template <class TOutputPixel, class TProbabilityPixel>
+AtlasRegistrationMethod<TOutputPixel, TProbabilityPixel>
+::AtlasRegistrationMethod() : m_WarpGrid(3, 0),
+ m_UseNonLinearInterpolation(true),
+ m_DoneRegistration(false),
+ m_Modified(false),
+ m_AtlasLinearTransformChoice("Affine"),
+ m_ImageLinearTransformChoice("Rigid"),
+ m_DebugLevel(0)
+{
+ m_InputImageTissueRegion = NULL;
+ m_InputSpatialObjectTissueRegion = NULL;
+}
+
+template <class TOutputPixel, class TProbabilityPixel>
+AtlasRegistrationMethod<TOutputPixel, TProbabilityPixel>
+::~AtlasRegistrationMethod()
+{