Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

ENH: Add DWIConvert to BRAINSStandAlone

  • Loading branch information...
commit 321c21eb2eee4631401f2723fbe458f9984441f2 1 parent 6ad5552
@Chaircrusher Chaircrusher authored
Showing with 6,607 additions and 84 deletions.
  1. +1 −6 BRAINSTools.cmake
  2. +1 −0  DWIConvert/CMakeLists.txt
  3. +43 −0 DWIConvert/DWIConvert/CMakeLists.txt
  4. +857 −0 DWIConvert/DWIConvert/DCMTKFileReader.h
  5. +2,180 −0 DWIConvert/DWIConvert/DWIConvert.cxx
  6. +137 −0 DWIConvert/DWIConvert/DWIConvert.xml
  7. +371 −0 DWIConvert/DWIConvert/DWIConvertUtils.h
  8. +307 −0 DWIConvert/DWIConvert/ExtendedTesting/CMakeLists.txt
  9. +269 −0 DWIConvert/DWIConvert/ExtendedTesting/DWICompare.cxx
  10. +29 −0 DWIConvert/DWIConvert/ExtendedTesting/DWICompare.xml
  11. +22 −0 DWIConvert/DWIConvert/ExtendedTesting/DWIConvertTest.cxx
  12. +307 −0 DWIConvert/DWIConvert/ExtendedTesting/DWISimpleCompare.cxx
  13. +36 −0 DWIConvert/DWIConvert/ExtendedTesting/DWISimpleCompare.xml
  14. +101 −0 DWIConvert/DWIConvert/ExtendedTesting/DicomToNrrdDWICompareTest.cmake
  15. +78 −0 DWIConvert/DWIConvert/ExtendedTesting/FSLTextFileCompare.cxx
  16. +42 −0 DWIConvert/DWIConvert/ExtendedTesting/FSLToNrrdTest.cmake
  17. +155 −0 DWIConvert/DWIConvert/ExtendedTesting/LOCAL_itkDifferenceImageFilter.h
  18. +265 −0 DWIConvert/DWIConvert/ExtendedTesting/LOCAL_itkDifferenceImageFilter.hxx
  19. +244 −0 DWIConvert/DWIConvert/ExtendedTesting/MIDAS.cmake
  20. +59 −0 DWIConvert/DWIConvert/ExtendedTesting/NrrdToFSLTest.cmake
  21. +114 −0 DWIConvert/DWIConvert/ExtendedTesting/NrrdToNIfTI.cxx
  22. +516 −0 DWIConvert/DWIConvert/ExtendedTesting/itkTestMain.h
  23. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDx.bval.md5
  24. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDx.bvec.md5
  25. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDx.nii.gz.md5
  26. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDx.nrrd.md5
  27. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDx.tar.gz.md5
  28. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDxBigEndian.tar.gz.md5
  29. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDxt.nrrd.md5
  30. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/GeSignaHDxt.tar.gz.md5
  31. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva1.nrrd.md5
  32. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva1.tar.gz.md5
  33. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva2.nrrd.md5
  34. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva2.tar.gz.md5
  35. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva3.nrrd.md5
  36. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva3.tar.gz.md5
  37. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva4.md5
  38. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva4.nrrd.md5
  39. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva6.md5
  40. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva6.nrrd.md5
  41. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva7.md5
  42. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchieva7.nrrd.md5
  43. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/PhilipsAchievaBigEndian1.tar.gz.md5
  44. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrio-Syngo2004A-1.md5
  45. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrio-Syngo2004A-1.nrrd.md5
  46. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrio-Syngo2004A-2.md5
  47. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrio-Syngo2004A-2.nrrd.md5
  48. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTim1.md5
  49. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTim1.nrrd.md5
  50. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTim2.md5
  51. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTim2.nrrd.md5
  52. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTim3.md5
  53. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTim3.nrrd.md5
  54. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensTrioTimBigEndian1.md5
  55. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensVerio.md5
  56. +1 −0  DWIConvert/DWIConvert/ExtendedTesting/keys/SiemensVerio.nrrd.md5
  57. +154 −0 DWIConvert/DWIConvert/FSLToNrrd.cxx
  58. +39 −0 DWIConvert/DWIConvert/FindLIBICONV.cmake
  59. +116 −0 DWIConvert/DWIConvert/NrrdToFSL.cxx
  60. +20 −0 DWIConvert/DWIConvert/StringContains.h
  61. +41 −0 DWIConvert/DWIConvert/Testing/CMakeLists.txt
  62. +14 −0 DWIConvert/DWIConvert/Testing/DWIConvertTest.cxx
  63. +8 −0 DWIConvert/README.md
  64. +15 −0 DWIConvert/dcmtkPatchScript.cmake
  65. +16 −0 DWIConvert/libjpeg_configure_step.cmake.in
  66. +16 −0 DWIConvert/libjpeg_make_step.cmake.in
  67. +0 −7 SuperBuild.cmake
  68. +0 −71 SuperBuild/External_DCMTK.cmake
View
7 BRAINSTools.cmake
@@ -102,6 +102,7 @@ set(brains_modulenames
BRAINSConstellationDetector
BRAINSABC
ConvertBetweenFileFormats
+ DWIConvert
)
if(USE_DebugImageViewer)
@@ -119,10 +120,4 @@ foreach(modulename ${brains_modulenames})
endif()
endforeach()
-# Slicer builds DWIConvert itself
-if(NOT INTEGRATE_WITH_SLICER AND USE_DWIConvert)
- add_subdirectory(${DWIConvert_SOURCE_DIR}
- ${CMAKE_CURRENT_BINARY_DIR}/DWIConvert)
-endif()
-
ExternalData_Add_Target( ${PROJECT_NAME}FetchData ) # Name of data management target
View
1  DWIConvert/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(DWIConvert)
View
43 DWIConvert/DWIConvert/CMakeLists.txt
@@ -0,0 +1,43 @@
+project (DWIConvert)
+
+#-----------------------------------------------------------------------------
+enable_testing()
+
+find_package(ITK 4 REQUIRED)
+
+set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR} ${CMAKE_MODULE_PATH})
+
+# SlicerExecutionModel
+find_package(SlicerExecutionModel REQUIRED GenerateCLP)
+include(${GenerateCLP_USE_FILE})
+
+if(${DCMTK_WITH_ICONV})
+find_package(LIBICONV REQUIRED)
+endif()
+
+set(DWIConvertSupportLib_SRCS
+ FSLToNrrd.cxx
+ NrrdToFSL.cxx
+ )
+
+set(DWIConvertTest_SRC)
+foreach(f ${DWIConvertSupportLib_SRCS})
+ list(APPEND DWIConvertTest_SRC
+ ${CMAKE_CURRENT_LIST_DIR}/${f})
+endforeach()
+
+add_library(DWIConvertSupportLib STATIC ${DWIConvertSupportLib_SRCS})
+target_link_libraries(DWIConvertSupportLib ${ITK_LIBRARIES} )
+
+set(MODULE_NAME DWIConvert)
+#-----------------------------------------------------------------------------
+SEMMacroBuildCLI(
+ NAME ${MODULE_NAME}
+ LOGO_HEADER ${Slicer_SOURCE_DIR}/Resources/NAMICLogo.h
+ TARGET_LIBRARIES DWIConvertSupportLib
+)
+
+# several files needed down in ExtenededTesting
+if(BUILD_TESTING)
+ add_subdirectory(ExtendedTesting)
+endif()
View
857 DWIConvert/DWIConvert/DCMTKFileReader.h
@@ -0,0 +1,857 @@
+#ifndef DCMTKFileReader_h
+#define DCMTKFileReader_h
+#undef HAVE_SSTREAM // 'twould be nice if people coded without using
+ // incredibly generic macro names
+#include "osconfig.h" // make sure OS specific configuration is included first
+
+#define INCLUDE_CSTDIO
+#define INCLUDE_CSTRING
+#include "ofstdinc.h"
+#include "dcvrds.h"
+#include "dcdict.h" // For DcmDataDictionary
+#include "dctk.h" /* for various dcmdata headers */
+#include "cmdlnarg.h" /* for prepareCmdLineArgs */
+#include "dcuid.h" /* for dcmtk version name */
+#include "dcrledrg.h" /* for DcmRLEDecoderRegistration */
+
+#include "dcmimage.h" /* for DicomImage */
+#include "digsdfn.h" /* for DiGSDFunction */
+#include "diciefn.h" /* for DiCIELABFunction */
+
+#include "ofconapp.h" /* for OFConsoleApplication */
+#include "ofcmdln.h" /* for OFCommandLine */
+
+#include "diregist.h" /* include to support color images */
+#include "ofstd.h" /* for OFStandard */
+#include "itkByteSwapper.h"
+#include "itkIntTypes.h"
+
+#define DCMTKException(body) \
+ { \
+ if(throwException) \
+ { \
+ itkGenericExceptionMacro(body); \
+ } \
+ else \
+ { \
+ return EXIT_FAILURE; \
+ } \
+ }
+// Don't print error messages if you're not throwing
+// an exception
+// std::cerr body;
+
+class DCMTKSequence
+{
+public:
+ DCMTKSequence() : m_DcmSequenceOfItems(0) {}
+ void SetDcmSequenceOfItems(DcmSequenceOfItems *seq)
+ {
+ this->m_DcmSequenceOfItems = seq;
+ }
+ int card() { return this->m_DcmSequenceOfItems->card(); }
+ int GetSequence(unsigned long index,
+ DCMTKSequence &target,bool throwException = true)
+ {
+ DcmItem *item = this->m_DcmSequenceOfItems->getItem(index);
+ DcmSequenceOfItems *sequence =
+ dynamic_cast<DcmSequenceOfItems *>(item);
+ if(sequence == 0)
+ {
+ DCMTKException(<< "Can't find DCMTKSequence at index " << index);
+ }
+ target.SetDcmSequenceOfItems(sequence);
+ return EXIT_SUCCESS;
+ }
+ int GetStack(unsigned short group,
+ unsigned short element,
+ DcmStack &resultStack, bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ if(this->m_DcmSequenceOfItems->search(tagkey,resultStack) != EC_Normal)
+ {
+ DCMTKException(<< "Can't find tag " << std::hex << group << " "
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+
+ int GetElementCS(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ DcmStack resultStack;
+ this->GetStack(group,element,resultStack);
+ DcmCodeString *codeStringElement = dynamic_cast<DcmCodeString *>(resultStack.top());
+ if(codeStringElement == 0)
+ {
+ DCMTKException(<< "Can't get CodeString Element at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ OFString ofString;
+ if(codeStringElement->getOFStringArray(ofString) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get OFString Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ target = "";
+ for(unsigned j = 0; j < ofString.length(); ++j)
+ {
+ target += ofString[j];
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementFD(unsigned short group,
+ unsigned short element,
+ double * &target,
+ bool throwException = true)
+ {
+ DcmStack resultStack;
+ this->GetStack(group,element,resultStack);
+ DcmFloatingPointDouble *fdItem = dynamic_cast<DcmFloatingPointDouble *>(resultStack.top());
+ if(fdItem == 0)
+ {
+ DCMTKException(<< "Can't get CodeString Element at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ if(fdItem->getFloat64Array(target) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get floatarray Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementFD(unsigned short group,
+ unsigned short element,
+ double &target,
+ bool throwException = true)
+ {
+ double *array;
+ this->GetElementFD(group,element,array,throwException);
+ target = array[0];
+ return EXIT_SUCCESS;
+ }
+ int GetElementDS(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ DcmStack resultStack;
+ this->GetStack(group,element,resultStack);
+ DcmDecimalString *decimalStringElement = dynamic_cast<DcmDecimalString *>(resultStack.top());
+ if(decimalStringElement == 0)
+ {
+ DCMTKException(<< "Can't get DecimalString Element at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ // check for # of expected numbers in DS
+ std::stringstream ss;
+ ss << target;
+ OFString arity(ss.str().c_str());
+ if(decimalStringElement->checkValue(arity) != EC_Normal)
+ {
+ DCMTKException(<< "Value doesn't have proper number of elements");
+ }
+ OFString ofString;
+ if(decimalStringElement->getOFStringArray(ofString) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get DecimalString Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ target = "";
+ for(unsigned j = 0; j < ofString.length(); ++j)
+ {
+ target += ofString[j];
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementSQ(unsigned short group,
+ unsigned short element,
+ DCMTKSequence &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmStack resultStack;
+ this->GetStack(group,element,resultStack);
+
+ DcmSequenceOfItems *seqElement = dynamic_cast<DcmSequenceOfItems *>(resultStack.top());
+ if(seqElement == 0)
+ {
+ DCMTKException(<< "Can't get at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ target.SetDcmSequenceOfItems(seqElement);
+ return EXIT_SUCCESS;
+ }
+private:
+ DcmSequenceOfItems *m_DcmSequenceOfItems;
+};
+
+class DCMTKFileReader
+{
+public:
+ typedef DCMTKFileReader Self;
+
+ DCMTKFileReader() : m_DFile(0),
+ m_Dataset(0),
+ m_Xfer(EXS_Unknown),
+ m_FrameCount(0),
+ m_FileNumber(-1L)
+ {
+ }
+ ~DCMTKFileReader()
+ {
+
+ delete m_DFile;
+ }
+ void SetFileName(const std::string &fileName)
+ {
+ this->m_FileName = fileName;
+ }
+ const std::string &GetFileName() const
+ {
+ return this->m_FileName;
+ }
+ void LoadFile()
+ {
+ if(this->m_FileName == "")
+ {
+ itkGenericExceptionMacro(<< "No filename given" );
+ }
+ if(this->m_DFile != 0)
+ {
+ delete this->m_DFile;
+ }
+ this->m_DFile = new DcmFileFormat();
+ OFCondition cond = this->m_DFile->loadFile(this->m_FileName.c_str());
+ if(cond.bad())
+ {
+ itkGenericExceptionMacro(<< cond.text() << ": reading file " << this->m_FileName);
+ }
+ this->m_Dataset = this->m_DFile->getDataset();
+ this->m_Xfer = this->m_Dataset->getOriginalXfer();
+ if(this->m_Dataset->findAndGetSint32(DCM_NumberOfFrames,this->m_FrameCount).bad())
+ {
+ this->m_FrameCount = 1;
+ }
+ int fnum;
+ this->GetElementIS(0x0020,0x0013,fnum);
+ this->m_FileNumber = fnum;
+ }
+ int GetElementLO(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmLongString *loItem = dynamic_cast<DcmLongString *>(el);
+ if(loItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ OFString ofString;
+ if(loItem->getOFStringArray(ofString) != EC_Normal)
+ {
+ DCMTKException(<< "Cant get string from element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ target = "";
+ for(unsigned i = 0; i < ofString.size(); i++)
+ {
+ target += ofString[i];
+ }
+ return EXIT_SUCCESS;
+ }
+
+ int GetElementLO(unsigned short group,
+ unsigned short element,
+ std::vector<std::string> &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmLongString *loItem = dynamic_cast<DcmLongString *>(el);
+ if(loItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ target.clear();
+ OFString ofString;
+ for(unsigned long i = 0; loItem->getOFString(ofString,i) == EC_Normal; ++i)
+ {
+ std::string targetStr = "";
+ for(unsigned j = 0; j < ofString.size(); j++)
+ {
+ targetStr += ofString[j];
+ }
+ target.push_back(targetStr);
+ }
+ return EXIT_SUCCESS;
+ }
+
+ /** Get an array of data values, as contained in a DICOM
+ * DecimalString Item
+ */
+ template <typename TType>
+ int GetElementDS(unsigned short group,
+ unsigned short element,
+ unsigned short count,
+ TType *target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmDecimalString *dsItem = dynamic_cast<DcmDecimalString *>(el);
+ if(dsItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ OFVector<Float64> doubleVals;
+ if(dsItem->getFloat64Vector(doubleVals) != EC_Normal)
+ {
+ DCMTKException(<< "Cant extract Array from DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(doubleVals.size() != count)
+ {
+ DCMTKException(<< "DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << " expected "
+ << count << "items, but found "
+ << doubleVals.size() << std::dec);
+
+ }
+ for(unsigned i = 0; i < count; i++)
+ {
+ target[i] = static_cast<TType>(doubleVals[i]);
+ }
+ return EXIT_SUCCESS;
+ }
+
+ template <typename TType>
+ int GetElementDSorOB(unsigned short group,
+ unsigned short element,
+ TType &target,
+ bool throwException = true)
+ {
+ if(this->GetElementDS<TType>(group,element,1,&target,false) == EXIT_SUCCESS)
+ {
+ return EXIT_SUCCESS;
+ }
+ std::string val;
+ if(this->GetElementOB(group,element,val) != EXIT_SUCCESS)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ const char *data = val.c_str();
+ const TType *fptr = reinterpret_cast<const TType *>(data);
+ target = *fptr;
+ switch(this->GetTransferSyntax())
+ {
+ case EXS_LittleEndianImplicit:
+ case EXS_LittleEndianExplicit:
+ itk::ByteSwapper<TType>::SwapFromSystemToLittleEndian(&target);
+ break;
+ case EXS_BigEndianImplicit:
+ case EXS_BigEndianExplicit:
+ itk::ByteSwapper<TType>::SwapFromSystemToBigEndian(&target);
+ break;
+ default:
+ break;
+ }
+ return EXIT_SUCCESS;
+
+ }
+ /** Get a DecimalString Item as a single string
+ */
+ int GetElementDS(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmDecimalString *dsItem = dynamic_cast<DcmDecimalString *>(el);
+ if(dsItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ OFString ofString;
+ if(dsItem->getOFStringArray(ofString) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get DecimalString Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ target = "";
+ for(unsigned j = 0; j < ofString.length(); ++j)
+ {
+ target += ofString[j];
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementFD(unsigned short group,
+ unsigned short element,
+ double &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmFloatingPointDouble *fdItem = dynamic_cast<DcmFloatingPointDouble *>(el);
+ if(fdItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(fdItem->getFloat64(target) != EC_Normal)
+ {
+ DCMTKException(<< "Cant extract Array from DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementFD(unsigned short group,
+ unsigned short element,
+ double * &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmFloatingPointDouble *fdItem = dynamic_cast<DcmFloatingPointDouble *>(el);
+ if(fdItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(fdItem->getFloat64Array(target) != EC_Normal)
+ {
+ DCMTKException(<< "Cant extract Array from DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementFL(unsigned short group,
+ unsigned short element,
+ float &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmFloatingPointSingle *flItem = dynamic_cast<DcmFloatingPointSingle *>(el);
+ if(flItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(flItem->getFloat32(target) != EC_Normal)
+ {
+ DCMTKException(<< "Cant extract Array from DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementFLorOB(unsigned short group,
+ unsigned short element,
+ float &target,
+ bool throwException = true)
+ {
+ if(this->GetElementFL(group,element,target,false) == EXIT_SUCCESS)
+ {
+ return EXIT_SUCCESS;
+ }
+ std::string val;
+ if(this->GetElementOB(group,element,val) != EXIT_SUCCESS)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ const char *data = val.c_str();
+ const float *fptr = reinterpret_cast<const float *>(data);
+ target = *fptr;
+ switch(this->GetTransferSyntax())
+ {
+ case EXS_LittleEndianImplicit:
+ case EXS_LittleEndianExplicit:
+ itk::ByteSwapper<float>::SwapFromSystemToLittleEndian(&target);
+ break;
+ case EXS_BigEndianImplicit:
+ case EXS_BigEndianExplicit:
+ itk::ByteSwapper<float>::SwapFromSystemToBigEndian(&target);
+ break;
+ default:
+ break;
+ }
+ return EXIT_SUCCESS;
+ }
+
+ int GetElementUS(unsigned short group,
+ unsigned short element,
+ unsigned short &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmUnsignedShort *usItem = dynamic_cast<DcmUnsignedShort *>(el);
+ if(usItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(usItem->getUint16(target) != EC_Normal)
+ {
+ DCMTKException(<< "Cant extract Array from DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+ int GetElementUS(unsigned short group,
+ unsigned short element,
+ unsigned short *&target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmUnsignedShort *usItem = dynamic_cast<DcmUnsignedShort *>(el);
+ if(usItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(usItem->getUint16Array(target) != EC_Normal)
+ {
+ DCMTKException(<< "Cant extract Array from DecimalString " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+ /** Get a DecimalString Item as a single string
+ */
+ int GetElementCS(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmCodeString *csItem = dynamic_cast<DcmCodeString *>(el);
+ if(csItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ OFString ofString;
+ if(csItem->getOFStringArray(ofString) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get DecimalString Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ target = "";
+ for(unsigned j = 0; j < ofString.length(); ++j)
+ {
+ target += ofString[j];
+ }
+ return EXIT_SUCCESS;
+ }
+
+ /** get an IS (Integer String Item
+ */
+ int GetElementIS(unsigned short group,
+ unsigned short element,
+ ::itk::int32_t &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmIntegerString *isItem = dynamic_cast<DcmIntegerString *>(el);
+ if(isItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ if(isItem->getSint32(target) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get DecimalString Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ return EXIT_SUCCESS;
+ }
+
+ int GetElementISorOB(unsigned short group,
+ unsigned short element,
+ ::itk::int32_t &target,
+ bool throwException = true)
+ {
+ if(this->GetElementIS(group,element,target,false) == EXIT_SUCCESS)
+ {
+ return EXIT_SUCCESS;
+ }
+ std::string val;
+ if(this->GetElementOB(group,element,val,throwException) != EXIT_SUCCESS)
+ {
+ return EXIT_FAILURE;
+ }
+ const char *data = val.c_str();
+ const int *iptr = reinterpret_cast<const int *>(data);
+ target = *iptr;
+ switch(this->GetTransferSyntax())
+ {
+ case EXS_LittleEndianImplicit:
+ case EXS_LittleEndianExplicit:
+ itk::ByteSwapper<int>::SwapFromSystemToLittleEndian(&target);
+ break;
+ case EXS_BigEndianImplicit:
+ case EXS_BigEndianExplicit:
+ itk::ByteSwapper<int>::SwapFromSystemToBigEndian(&target);
+ break;
+ default: // no idea what to do
+ break;
+ }
+
+ return EXIT_SUCCESS;
+ }
+
+ int GetElementCSorOB(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ if(this->GetElementCS(group,element,target,false) == EXIT_SUCCESS)
+ {
+ return EXIT_SUCCESS;
+ }
+ std::string val;
+ if(this->GetElementOB(group,element,target,throwException) != EXIT_SUCCESS)
+ {
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+ }
+
+ /** get an OB OtherByte Item
+ */
+ int GetElementOB(unsigned short group,
+ unsigned short element,
+ std::string &target,
+ bool throwException = true)
+ {
+ DcmTagKey tagkey(group,element);
+ DcmElement *el;
+ if(this->m_Dataset->findAndGetElement(tagkey,el) != EC_Normal)
+ {
+ DCMTKException(<< "Cant find tag " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ DcmOtherByteOtherWord *obItem = dynamic_cast<DcmOtherByteOtherWord *>(el);
+ if(obItem == 0)
+ {
+ DCMTKException(<< "Cant find DecimalString element " << std::hex
+ << group << " " << std::hex
+ << element << std::dec);
+ }
+ OFString ofString;
+ if(obItem->getOFStringArray(ofString) != EC_Normal)
+ {
+ DCMTKException(<< "Can't get OFString Value at tag "
+ << std::hex << group << " "
+ << element << std::dec);
+ }
+ target = Self::ConvertFromOB(ofString);
+ return EXIT_SUCCESS;
+ }
+
+ int GetElementSQ(unsigned short group,
+ unsigned short entry,
+ DCMTKSequence &sequence,
+ bool throwException = true)
+ {
+ DcmSequenceOfItems *seq;
+ DcmTagKey tagKey(group,entry);
+
+ if(this->m_Dataset->findAndGetSequence(tagKey,seq) != EC_Normal)
+ {
+ DCMTKException(<< "Can't find sequence "
+ << std::hex << group << " "
+ << std::hex << entry)
+ }
+ sequence.SetDcmSequenceOfItems(seq);
+ return EXIT_SUCCESS;
+ }
+ int GetFrameCount() { return this->m_FrameCount; }
+
+ E_TransferSyntax GetTransferSyntax() { return m_Xfer; }
+
+ long GetFileNumber() const
+ {
+ return m_FileNumber;
+ }
+ static void
+ AddDictEntry(DcmDictEntry *entry)
+ {
+ DcmDataDictionary &dict = dcmDataDict.wrlock();
+ dict.addEntry(entry);
+ dcmDataDict.unlock();
+ }
+
+private:
+ static unsigned ascii2hex(char c)
+ {
+ switch(c)
+ {
+ case '0': return 0;
+ case '1': return 1;
+ case '2': return 2;
+ case '3': return 3;
+ case '4': return 4;
+ case '5': return 5;
+ case '6': return 6;
+ case '7': return 7;
+ case '8': return 8;
+ case '9': return 9;
+ case 'a':
+ case 'A': return 10;
+ case 'b':
+ case 'B': return 11;
+ case 'c':
+ case 'C': return 12;
+ case 'd':
+ case 'D': return 13;
+ case 'e':
+ case 'E': return 14;
+ case 'f':
+ case 'F': return 15;
+ }
+ return 255; // should never happen
+ }
+ static std::string ConvertFromOB(OFString &toConvert)
+ {
+ // string format is nn\nn\nn...
+ std::string rval;
+ for(size_t pos = 0; pos < toConvert.size(); pos += 3)
+ {
+ unsigned char convert[2];
+ convert[0] = Self::ascii2hex(toConvert[pos]);
+ convert[1] = Self::ascii2hex(toConvert[pos+1]);
+ unsigned char conv = convert[0] << 4;
+ conv += convert[1];
+ rval.push_back(static_cast<unsigned char>(conv));
+ }
+ return rval;
+ }
+
+ std::string m_FileName;
+ DcmFileFormat* m_DFile;
+ DcmDataset * m_Dataset;
+ E_TransferSyntax m_Xfer;
+ Sint32 m_FrameCount;
+ long m_FileNumber;
+};
+
+bool CompareDCMTKFileReaders(DCMTKFileReader *a, DCMTKFileReader *b)
+{
+ return a->GetFileNumber() < b->GetFileNumber();
+}
+
+#endif // DCMTKFileReader_h
View
2,180 DWIConvert/DWIConvert/DWIConvert.cxx
@@ -0,0 +1,2180 @@
+/*=========================================================================
+Computing (NAMIC), funded by the National Institutes of Health
+through the NIH Roadmap for Medical Research, Grant U54 EB005149.
+
+See License.txt or http://www.slicer.org/copyright/copyright.txt for details.
+
+ ***
+ This program converts Diffusion weighted MR images in Dicom format into
+ NRRD format.
+
+Assumptions:
+
+1) Uses left-posterior-superior (Dicom default) as default space for philips and siemens.
+This is the default space for NRRD header.
+2) For GE data, Dicom data are arranged in volume interleaving order.
+3) For Siemens data, images are arranged in mosaic form.
+4) For oblique collected Philips data, the measurement frame for the
+gradient directions is the same as the ImageOrientationPatient
+
+Reference materials:
+DICOM Data Dictionary: http://medical.nema.org/Dicom/2011/11_06pu.pdf
+=========================================================================*/
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <algorithm>
+#include <ctype.h>
+#include "DWIConvertCLP.h"
+
+#include "itkMacro.h"
+#include "itkIntTypes.h"
+#include "itkDCMTKSeriesFileNames.h"
+#undef HAVE_SSTREAM
+#include "itkDCMTKImageIO.h"
+#include "itkRawImageIO.h"
+#include "itkImage.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkImageSeriesReader.h"
+#include "itksys/Directory.hxx"
+#include "itksys/SystemTools.hxx"
+#include "itksys/Base64.h"
+#undef HAVE_SSTREAM
+#include "itkDCMTKFileReader.h"
+#include "djdecode.h"
+#include "StringContains.h"
+#include "DWIConvertUtils.h"
+
+unsigned int ConvertFromCharPtr(const char *s)
+{
+ unsigned int rval = 0;
+ // assume little-endian
+ for(unsigned i = 0; i < sizeof(unsigned int); ++i)
+ {
+ rval += ((unsigned int)s[i]) << (i * 8);
+ }
+// this makes no sense, according to what I've read, but apparently,
+// the uint32 numbers in the CSA header are little-endian even in
+// files with BigEndian transfer syntax.#if 0
+#if 0
+ switch(xferSyntax)
+ {
+ case EXS_LittleEndianImplicit:
+ case EXS_LittleEndianExplicit:
+ itk::ByteSwapper<unsigned int>::SwapFromSystemToLittleEndian(&rval);
+ break;
+ case EXS_BigEndianImplicit:
+ case EXS_BigEndianExplicit:
+ itk::ByteSwapper<unsigned int>::SwapFromSystemToBigEndian(&rval);
+ break;
+ default:
+ break;
+ }
+#endif
+ return rval;
+}
+/** pull data out of Siemens scans.
+ *
+ * Siemens sticks most of the DTI information into a single
+ * OB-format entry. This is actually rigidly structured, but
+ * this function depends on the needed data living at fixed offset
+ * from the beginning of the name of each tag, and ignores the
+ * internal structure documented in the Siemens Dicom Compliance
+ * document.
+ */
+unsigned int
+ExtractSiemensDiffusionInformation(const std::string tagString,
+ const std::string nameString,
+ std::vector<double>& valueArray)
+{
+ ::size_t atPosition = tagString.find( nameString );
+ if(atPosition == std::string::npos)
+ {
+ return 0;
+ }
+ while( true ) // skip nameString inside a quotation
+ {
+ std::string nextChar = tagString.substr( atPosition+nameString.size(), 1 );
+
+ if (nextChar.c_str()[0] == 0 )
+ {
+ break;
+ }
+ else
+ {
+ atPosition = tagString.find( nameString, atPosition+2 );
+ }
+ }
+
+ if ( atPosition == std::string::npos)
+ {
+ return 0;
+ }
+ std::string infoAsString = tagString.substr( atPosition, tagString.size()-atPosition+1 );
+ const char * infoAsCharPtr = infoAsString.c_str();
+
+ unsigned int vm = ConvertFromCharPtr(infoAsCharPtr+64);
+ {
+ std::string vr = infoAsString.substr( 68, 2 );
+ int syngodt = ConvertFromCharPtr(infoAsCharPtr+72);
+ int nItems = ConvertFromCharPtr(infoAsCharPtr+76);
+ int localDummy = ConvertFromCharPtr(infoAsCharPtr+80);
+
+ //std::cout << "\tName String: " << nameString << std::endl;
+ //std::cout << "\tVR: " << vr << std::endl;
+ //std::cout << "\tVM: " << vm << std::endl;
+ //std::cout << "Local String: " << infoAsString.substr(0,80) << std::endl;
+
+ /* This hack is required for some Siemens VB15 Data */
+ if ( ( nameString == "DiffusionGradientDirection" ) && (vr != "FD") )
+ {
+ bool loop = true;
+ while ( loop )
+ {
+ atPosition = tagString.find( nameString, atPosition+26 );
+ if ( atPosition == std::string::npos)
+ {
+ //std::cout << "\tFailed to find DiffusionGradientDirection Tag - returning" << vm << std::endl;
+ return 0;
+ }
+ infoAsString = tagString.substr( atPosition, tagString.size()-atPosition+1 );
+ infoAsCharPtr = infoAsString.c_str();
+ //std::cout << "\tOffset to new position" << std::endl;
+ //std::cout << "\tNew Local String: " << infoAsString.substr(0,80) << std::endl;
+ vm = ConvertFromCharPtr(infoAsCharPtr+64);
+ vr = infoAsString.substr( 68, 2 );
+ if (vr == "FD") loop = false;
+ syngodt = ConvertFromCharPtr(infoAsCharPtr+72);
+ nItems = ConvertFromCharPtr(infoAsCharPtr+76);
+ localDummy = ConvertFromCharPtr(infoAsCharPtr+80);
+ //std::cout << "\tVR: " << vr << std::endl;
+ //std::cout << "\tVM: " << vm << std::endl;
+ }
+ }
+ else
+ {
+ //std::cout << "\tUsing initial position" << std::endl;
+ }
+ //std::cout << "\tArray Length: " << vm << std::endl;
+ }
+
+ unsigned int offset = 84;
+ for (unsigned int k = 0; k < vm; ++k)
+ {
+ const int itemLength = ConvertFromCharPtr(infoAsCharPtr+offset+4);
+ const int strideSize = static_cast<int> (ceil(static_cast<double>(itemLength)/4) * 4);
+ if( infoAsString.length() < offset + 16 + itemLength )
+ {
+ // data not available or incomplete
+ return 0;
+ }
+ const std::string valueString = infoAsString.substr( offset+16, itemLength );
+ valueArray.push_back( atof(valueString.c_str()) );
+ offset += 16+strideSize;
+ }
+ return vm;
+}
+
+/**
+ * Add private tags to the Dicom Dictionary
+ */
+void
+AddFlagsToDictionary()
+{
+ // these have to be dynamically allocated because otherwise there's
+ // a malloc error after main exits.
+
+ // relevant GE tags
+ DcmDictEntry *GEDictBValue = new DcmDictEntry(0x0043, 0x1039, DcmVR(EVR_IS),
+ "B Value of diffusion weighting", 1, 1, 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *GEDictXGradient = new DcmDictEntry(0x0019, 0x10bb, DcmVR(EVR_DS),
+ "X component of gradient direction", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *GEDictYGradient = new DcmDictEntry(0x0019, 0x10bc, DcmVR(EVR_DS),
+ "Y component of gradient direction", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *GEDictZGradient = new DcmDictEntry(0x0019, 0x10bd, DcmVR(EVR_DS),
+ "Z component of gradient direction", 1, 1 , 0,true,
+ "dicomtonrrd");
+
+ // relevant Siemens private tags
+ DcmDictEntry *SiemensMosiacParameters = new DcmDictEntry(0x0051, 0x100b, DcmVR(EVR_IS),
+ "Mosiac Matrix Size", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *SiemensDictNMosiac = new DcmDictEntry(0x0019, 0x100a, DcmVR(EVR_US),
+ "Number of Images In Mosaic", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *SiemensDictBValue = new DcmDictEntry(0x0019, 0x100c, DcmVR(EVR_IS),
+ "B Value of diffusion weighting", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *SiemensDictDiffusionDirection = new DcmDictEntry(0x0019, 0x100e, DcmVR(EVR_FD),
+ "Diffusion Gradient Direction", 3, 3 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *SiemensDictDiffusionMatrix = new DcmDictEntry(0x0019, 0x1027, DcmVR(EVR_FD),
+ "Diffusion Matrix", 6, 6 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *SiemensDictShadowInfo = new DcmDictEntry(0x0029, 0x1010, DcmVR(EVR_OB),
+ "Siemens DWI Info", 1, 1 , 0,true,
+ "dicomtonrrd");
+
+ // relevant Philips private tags
+ DcmDictEntry *PhilipsDictBValue = new DcmDictEntry(0x2001, 0x1003, DcmVR(EVR_FL),
+ "B Value of diffusion weighting", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *PhilipsDictDiffusionDirection = new DcmDictEntry(0x2001, 0x1004, DcmVR(EVR_CS),
+ "Diffusion Gradient Direction", 1, 1 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *PhilipsDictDiffusionDirectionRL = new DcmDictEntry(0x2005, 0x10b0, DcmVR(EVR_FL),
+ "Diffusion Direction R/L", 4, 4 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *PhilipsDictDiffusionDirectionAP = new DcmDictEntry(0x2005, 0x10b1, DcmVR(EVR_FL),
+ "Diffusion Direction A/P", 4, 4 , 0,true,
+ "dicomtonrrd");
+ DcmDictEntry *PhilipsDictDiffusionDirectionFH = new DcmDictEntry(0x2005, 0x10b2, DcmVR(EVR_FL),
+ "Diffusion Direction F/H", 4, 4 , 0,true,
+ "dicomtonrrd");
+
+ itk::DCMTKFileReader::AddDictEntry(GEDictBValue);
+ itk::DCMTKFileReader::AddDictEntry(GEDictXGradient);
+ itk::DCMTKFileReader::AddDictEntry(GEDictYGradient);
+ itk::DCMTKFileReader::AddDictEntry(GEDictZGradient);
+
+ // relevant Siemens private tags
+ itk::DCMTKFileReader::AddDictEntry(SiemensMosiacParameters);
+ itk::DCMTKFileReader::AddDictEntry(SiemensDictNMosiac);
+ itk::DCMTKFileReader::AddDictEntry(SiemensDictBValue);
+ itk::DCMTKFileReader::AddDictEntry(SiemensDictDiffusionDirection);
+ itk::DCMTKFileReader::AddDictEntry(SiemensDictDiffusionMatrix);
+ itk::DCMTKFileReader::AddDictEntry(SiemensDictShadowInfo);
+
+ // relevant Philips private tags
+ itk::DCMTKFileReader::AddDictEntry(PhilipsDictBValue);
+ itk::DCMTKFileReader::AddDictEntry(PhilipsDictDiffusionDirection);
+ itk::DCMTKFileReader::AddDictEntry(PhilipsDictDiffusionDirectionRL);
+ itk::DCMTKFileReader::AddDictEntry(PhilipsDictDiffusionDirectionAP);
+ itk::DCMTKFileReader::AddDictEntry(PhilipsDictDiffusionDirectionFH);
+
+}
+
+/** Free the headers for all dicom files.
+ * also calls the DJDecoder cleanup.
+ */
+void FreeHeaders(std::vector<itk::DCMTKFileReader *> &allHeaders)
+{
+ for(std::vector<itk::DCMTKFileReader *>::iterator it = allHeaders.begin();
+ it != allHeaders.end(); ++it)
+ {
+ delete (*it);
+ }
+}
+
+typedef short PixelValueType;
+typedef itk::Image< PixelValueType, 3 > VolumeType;
+
+int
+Write4DVolume( VolumeType::Pointer &img, int nVolumes, const std::string &fname )
+{
+ typedef itk::Image<PixelValueType,4> Volume4DType;
+
+ VolumeType::SizeType size3D(img->GetLargestPossibleRegion().GetSize());
+ VolumeType::DirectionType direction3D(img->GetDirection());
+ VolumeType::SpacingType spacing3D(img->GetSpacing());
+ VolumeType::PointType origin3D(img->GetOrigin());
+
+ Volume4DType::SizeType size4D;
+ size4D[0] = size3D[0];
+ size4D[1] = size3D[1];
+ size4D[2] = size3D[2] / nVolumes;
+ size4D[3] = nVolumes;
+
+ Volume4DType::DirectionType direction4D;
+ Volume4DType::SpacingType spacing4D;
+ Volume4DType::PointType origin4D;
+
+ for(unsigned i = 0; i < 3; ++i)
+ {
+ for(unsigned j = 0; j < 3; ++j)
+ {
+ direction4D[i][j] = direction3D[i][j];
+ }
+ direction4D[3][i] = 0.0;
+ direction4D[i][3] = 0.0;
+ spacing4D[i] = spacing3D[i];
+ origin4D[i] = origin3D[i];
+ }
+ direction4D[3][3] = 1.0;
+ spacing4D[3] = 1.0;
+ origin4D[3] = 0.0;
+
+
+ Volume4DType::Pointer img4D = Volume4DType::New();
+ img4D->SetRegions(size4D);
+ img4D->SetDirection(direction4D);
+ img4D->SetSpacing(spacing4D);
+ img4D->SetOrigin(origin4D);
+
+ img4D->Allocate();
+ memcpy(img4D->GetBufferPointer(),
+ img->GetBufferPointer(),
+ img->GetLargestPossibleRegion().GetNumberOfPixels() * sizeof(PixelValueType));
+#if 0
+ {
+ itk::ImageFileWriter< VolumeType >::Pointer writer = itk::ImageFileWriter< VolumeType >::New();
+ writer->SetFileName( "dwi3dconvert.nii.gz");
+ writer->SetInput( img );
+ writer->Update();
+ }
+#endif
+
+ itk::ImageFileWriter< Volume4DType >::Pointer imgWriter =
+ itk::ImageFileWriter< Volume4DType >::New();
+
+ imgWriter->SetInput( img4D );
+ imgWriter->SetFileName( fname.c_str() );
+ try
+ {
+ imgWriter->Update();
+ }
+ catch (itk::ExceptionObject &excp)
+ {
+ std::cerr << "Exception thrown while writing "
+ << fname << std::endl;
+ std::cerr << excp << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
+
+void
+DeInterleaveVolume(VolumeType::Pointer &volume,
+ size_t SlicesPerVolume,
+ size_t NSlices)
+{
+ size_t NVolumes = NSlices / SlicesPerVolume;
+
+ VolumeType::RegionType R = volume->GetLargestPossibleRegion();
+ R.SetSize(2,1);
+ std::vector<VolumeType::PixelType> v(NSlices);
+ std::vector<VolumeType::PixelType> w(NSlices);
+
+ itk::ImageRegionIteratorWithIndex<VolumeType> I( volume, R );
+ // permute the slices by extracting the 1D array of voxels for
+ // a particular {x,y} position, then re-ordering the voxels such
+ // that all the voxels for a particular volume are adjacent
+ for (I.GoToBegin(); !I.IsAtEnd(); ++I)
+ {
+ VolumeType::IndexType idx = I.GetIndex();
+
+ // extract all values in one "column"
+ for (unsigned int k = 0; k < NSlices; ++k)
+ {
+ idx[2] = k;
+ v[k] = volume->GetPixel( idx );
+ }
+
+ // permute
+ for (unsigned int k = 0; k < NVolumes; ++k)
+ {
+ for (unsigned int m = 0; m < SlicesPerVolume; ++m)
+ {
+ w[(k * SlicesPerVolume) + m] = v[ (m * NVolumes) + k];
+ }
+ }
+
+ // put things back in order
+ for (unsigned int k = 0; k < NSlices; ++k)
+ {
+ idx[2] = k;
+ volume->SetPixel( idx, w[k] );
+ }
+ }
+
+}
+
+int main(int argc, char *argv[])
+{
+ PARSE_ARGS;
+
+ if(conversionMode == "FSLToNrrd")
+ {
+ extern int FSLToNrrd(const std::string &inputVolume,
+ const std::string &outputVolume,
+ const std::string &inputBValues,
+ const std::string &inputBVectors);
+
+ return FSLToNrrd(inputVolume, outputVolume,
+ inputBValues, inputBVectors);
+ }
+ if(conversionMode == "NrrdToFSL")
+ {
+ extern int NrrdToFSL(const std::string &inputVolume,
+ const std::string &outputVolume,
+ const std::string &outputBValues,
+ const std::string &outputBVectors);
+ return NrrdToFSL(inputVolume,outputVolume,
+ outputBValues,outputBVectors);
+ }
+
+
+ typedef itk::ImageSeriesReader< VolumeType > ReaderType;
+ typedef itk::ImageFileReader< VolumeType > SingleFileReaderType;
+ typedef itk::DCMTKSeriesFileNames InputNamesGeneratorType;
+
+ AddFlagsToDictionary();
+ bool nrrdFormat(true);
+ //
+ // check for required parameters
+ if(inputDicomDirectory == "")
+ {
+ std::cerr << "Missing DICOM input directory path" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ if(outputVolume == "")
+ {
+ std::cerr << "Missing DICOM output volume name" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ std::string outputVolumeHeaderName(outputVolume);
+ if(outputVolume.find("/") == std::string::npos &&
+ outputVolume.find("\\") == std::string::npos)
+ {
+ if(outputVolumeHeaderName.size() != 0)
+ {
+ outputVolumeHeaderName = outputDirectory;
+ outputVolumeHeaderName += "/";
+ outputVolumeHeaderName += outputVolume;
+ }
+ }
+
+ // decide whether the output is a single file or
+ // header/raw pair
+ std::string outputVolumeDataName;
+ std::string outputFSLBValFilename;
+ std::string outputFSLBVecFilename;
+ if(conversionMode != "DicomToFSL")
+ {
+ const size_t extensionPos = outputVolumeHeaderName.find(".nhdr");
+ if(extensionPos != std::string::npos)
+ {
+ outputVolumeDataName = outputVolumeHeaderName.substr(0,extensionPos);
+ outputVolumeDataName += ".raw";
+ nrrdFormat = false;
+ }
+ }
+ else
+ {
+ // FSL output of gradients & BValues
+ std::string fslPrefix;
+ size_t extensionPos;
+ extensionPos = outputVolumeHeaderName.find(".nii.gz");
+ if(extensionPos == std::string::npos)
+ {
+ extensionPos = outputVolumeHeaderName.find(".nii");
+ if(extensionPos == std::string::npos)
+ {
+ std::cerr << "FSL Format output chosen, "
+ << "but output Volume not a recognized "
+ << "NIfTI filename " << outputVolumeHeaderName
+ << std::endl;
+ exit(1);
+ }
+ }
+ if(outputBValues == "")
+ {
+ outputFSLBValFilename = outputVolumeHeaderName.substr(0,extensionPos);
+ outputFSLBValFilename += ".bval";
+ }
+ else
+ {
+ outputFSLBValFilename = outputBValues;
+ }
+ if(outputBVectors == "")
+ {
+ outputFSLBVecFilename = outputVolumeHeaderName.substr(0,extensionPos);
+ outputFSLBVecFilename += ".bvec";
+ }
+ else
+ {
+ outputFSLBVecFilename = outputBVectors;
+ }
+ }
+
+ ReaderType::FileNamesContainer inputFileNames;
+ //
+ // get the names of all slices in the directory
+ InputNamesGeneratorType::Pointer inputNames = InputNamesGeneratorType::New();
+ if(itksys::SystemTools::FileIsDirectory(inputDicomDirectory.c_str()))
+ {
+ inputNames->SetUseSeriesDetails( true);
+ inputNames->SetLoadSequences( true );
+ inputNames->SetLoadPrivateTags( true );
+ inputNames->SetInputDirectory(inputDicomDirectory);
+ inputFileNames = inputNames->GetInputFileNames();
+ }
+ else if(itksys::SystemTools::FileExists(inputDicomDirectory.c_str()))
+ // or, if it isn't a directory, maybe it is
+ // a multi-frame DICOM file.
+ {
+ inputFileNames.push_back(inputDicomDirectory);
+ }
+
+ const size_t nFiles = inputFileNames.size();
+
+ if(nFiles < 1)
+ {
+ std::cerr << "Error: no DICOMfiles found in inputDirectory: " << inputDicomDirectory
+ << std::endl;
+ return EXIT_FAILURE;
+ }
+#if 0
+ // this duplicates what the itk::DCMTKSeriesFileNames does, so
+ // comment it out.
+ else if(inputFileNames.size() == 1)
+ {
+ // Use itksys::Directory to open up the given name as a directory.
+ inputFileNames.resize( 0 );
+ itksys::Directory directory;
+ directory.Load( itksys::SystemTools::CollapseFullPath(inputDicomDirectory.c_str()).c_str() );
+ typedef itk::DCMTKImageIO ImageIOType;
+
+ // for each patient directory
+ for ( unsigned int k = 0; k < directory.GetNumberOfFiles(); ++k)
+ {
+ std::string subdirectory( inputDicomDirectory.c_str() );
+ subdirectory = subdirectory + "/" + directory.GetFile(k);
+
+ const std::string sqDir( directory.GetFile(k) );
+ if (sqDir.length() == 1 && directory.GetFile(k)[0] == '.') // skip self
+ {
+ continue;
+ }
+ else if (sqDir.length() == 2 && sqDir.find( ".." ) != std::string::npos) // skip parent
+ {
+ continue;
+ }
+ else if (!itksys::SystemTools::FileIsDirectory( subdirectory.c_str() )) // load only files
+ {
+ if (itk::DCMTKFileReader::IsImageFile(subdirectory.c_str()))
+ {
+ inputFileNames.push_back( subdirectory );
+ }
+ }
+ }
+ }
+#endif
+
+ //////////////////////////////////////////////////
+ // load all files in the dicom series.
+ //////////////////////////////////////////////////
+
+ std::vector<itk::DCMTKFileReader *> allHeaders(inputFileNames.size());
+ int headerCount = 0;
+
+ for(unsigned i = 0; i < allHeaders.size(); ++i)
+ {
+ itk::DCMTKFileReader *curReader = new itk::DCMTKFileReader;
+ curReader->SetFileName(inputFileNames[i]);
+ try
+ {
+ curReader->LoadFile();
+ }
+ catch(...)
+ {
+ std::cerr << "Error reading slice" << inputFileNames[i] << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+ //
+ // check for pixel data.
+ if(!curReader->HasPixelData())
+ {
+ delete curReader;
+ }
+ else
+ {
+ allHeaders[headerCount] = curReader;
+ headerCount++;
+ }
+ }
+
+ if(headerCount == 0)
+ {
+ std::cerr << "No pixel data in series" << std::endl;
+ return EXIT_FAILURE;
+ }
+ //
+ // reorder the filename list
+ inputFileNames.resize( 0 );
+ for(unsigned i = 0; i < allHeaders.size(); ++i)
+ {
+ std::string fname = allHeaders[i]->GetFileName();
+ inputFileNames.push_back(fname);
+ }
+
+ //
+ // find vendor and modality
+ std::string vendor;
+ try
+ {
+ allHeaders[0]->GetElementLO(0x0008,0x0070,vendor);
+ strupper(vendor);
+ }
+ catch(itk::ExceptionObject &excp)
+ {
+ std::cerr << "Can't get vendor name from DICOM file" << excp << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+
+ std::string modality;
+ try
+ {
+ allHeaders[0]->GetElementCS(0x0008,0x0060,modality);
+ }
+ catch(itk::ExceptionObject &excp)
+ {
+ std::cerr << "Can't find modality in DICOM file" << excp << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+
+ //
+ // IF it's a PET or SPECT file, just write it out as a float image.
+ if(StringContains(modality,"PT") || StringContains(modality,"ST"))
+ {
+ typedef itk::Image<float, 3> USVolumeType;
+ itk::ImageSeriesReader<USVolumeType>::Pointer seriesReader =
+ itk::ImageSeriesReader<USVolumeType>::New();
+ seriesReader->SetFileNames( inputFileNames );
+
+ itk::ImageFileWriter<USVolumeType>::Pointer nrrdImageWriter =
+ itk::ImageFileWriter<USVolumeType>::New();
+
+ nrrdImageWriter->SetFileName( outputVolumeHeaderName );
+ nrrdImageWriter->SetInput( seriesReader->GetOutput() );
+ try
+ {
+ nrrdImageWriter->Update();
+ }
+ catch( itk::ExceptionObject & err )
+ {
+ std::cerr << "ExceptionObject caught !" << std::endl;
+ std::cerr << err << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+ FreeHeaders(allHeaders);
+ return EXIT_SUCCESS;
+ }
+
+ //
+ // any failure in extracting data is an error
+ // so encapsulate rest of program in exception block
+ try
+ {
+
+ bool SliceMosaic(false);
+ // decide if it is a mosaic
+ if(StringContains(vendor,"SIEMENS"))
+ {
+ std::string ImageType;
+ allHeaders[0]->GetElementCS(0x0008,0x0008, ImageType);
+ if(ImageType.find("MOSAIC") != std::string::npos)
+ {
+ SliceMosaic = true;
+ }
+ }
+ else if(!StringContains(vendor,"GE") && !StringContains(vendor,"PHILIPS"))
+ {
+ std::cerr << "Unrecognized scanner vendor |"
+ << vendor << "|" << std::endl;
+ }
+
+ //////////////////////////////////////////////////
+ // 1) Read the input series as an array of slices
+ unsigned int nSlice;
+ VolumeType::Pointer readerOutput;
+ itk::DCMTKImageIO::Pointer dcmtkIO = itk::DCMTKImageIO::New();
+ bool multiSliceVolume;
+ if(inputFileNames.size() > 1)
+ {
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetImageIO( dcmtkIO );
+ reader->SetFileNames( inputFileNames );
+ nSlice
+ = inputFileNames.size();
+ try
+ {
+ reader->Update();
+ }
+ catch (itk::ExceptionObject &excp)
+ {
+ std::cerr << "Exception thrown while reading the series" << std::endl;
+ std::cerr << excp << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+ readerOutput = reader->GetOutput();
+ multiSliceVolume = false;
+#if 0
+ itk::ImageFileWriter< VolumeType >::Pointer writer = itk::ImageFileWriter< VolumeType >::New();
+ writer->SetFileName( "dwiconvert.nrrd");
+ writer->SetInput( readerOutput );
+ writer->Update();
+#endif
+ }
+ else
+ {
+ SingleFileReaderType::Pointer reader =
+ SingleFileReaderType::New();
+ reader->SetImageIO( dcmtkIO );
+ reader->SetFileName( inputFileNames[0] );
+ nSlice = inputFileNames.size();
+ try
+ {
+ reader->Update();
+ }
+ catch (itk::ExceptionObject &excp)
+ {
+ std::cerr << "Exception thrown while reading the series" << std::endl;
+ std::cerr << excp << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+ readerOutput = reader->GetOutput();
+ multiSliceVolume = true;
+ }
+
+ // get image dims and resolution
+ unsigned short nRows, nCols;
+
+ allHeaders[0]->GetElementUS(0x0028,0x0010,nRows);
+ allHeaders[0]->GetElementUS(0x0028,0x0011,nCols);
+
+ double sliceSpacing;
+ double xRes, yRes;
+ {
+ double spacing[3];
+ allHeaders[0]->GetSpacing(spacing);
+ yRes = spacing[1];
+ xRes = spacing[0];
+ sliceSpacing = spacing[2];
+ }
+
+ itk::Vector<double,3> ImageOrigin;
+ {
+ double origin[3];
+ allHeaders[0]->GetOrigin(origin);
+ ImageOrigin[0] = origin[0];
+ ImageOrigin[1] = origin[1];
+ ImageOrigin[2] = origin[2];
+ }
+
+ unsigned int numberOfSlicesPerVolume;
+ std::map<std::string,int> sliceLocations;
+ if(!multiSliceVolume)
+ {
+ // Make a hash of the sliceLocations in order to get the correct
+ // count. This is more reliable since SliceLocation may not be available.
+ std::vector<int> sliceLocationIndicator;
+ std::vector<std::string> sliceLocationStrings;
+
+ sliceLocationIndicator.resize( nSlice );
+
+ for (unsigned int k = 0; k < nSlice; ++k)
+ {
+ std::string originString;
+
+ allHeaders[k]->GetElementDS(0x0020, 0x0032, originString );
+ sliceLocationStrings.push_back( originString );
+ sliceLocations[originString]++;
+ // std::cerr << inputFileNames[k] << " " << originString << std::endl;
+ }
+
+ for (unsigned int k = 0; k < nSlice; ++k)
+ {
+ std::map<std::string,int>::iterator it = sliceLocations.find( sliceLocationStrings[k] );
+ sliceLocationIndicator[k] = distance( sliceLocations.begin(), it );
+ }
+
+ numberOfSlicesPerVolume=sliceLocations.size();
+ std::cout << "=================== numberOfSlicesPerVolume:" << numberOfSlicesPerVolume << std::endl;
+
+ //
+ // if the numberOfSlicesPerVolume == 1, de-interleaving won't do
+ // anything so there's no point in doing it.
+ if ( nSlice >= 2 && numberOfSlicesPerVolume > 1)
+ {
+ if(sliceLocationIndicator[0] != sliceLocationIndicator[1])
+ {
+ std::cout << "Dicom images are ordered in a volume interleaving way." << std::endl;
+ }
+ else
+ {
+ std::cout << "Dicom images are ordered in a slice interleaving way." << std::endl;
+ // reorder slices into a volume interleaving manner
+ DeInterleaveVolume(readerOutput,numberOfSlicesPerVolume,nSlice);
+#if 0
+ itk::ImageFileWriter< VolumeType >::Pointer writer = itk::ImageFileWriter< VolumeType >::New();
+ writer->SetFileName( "deinterleave.nrrd");
+ writer->SetInput( readerOutput );
+ writer->Update();
+#endif
+ }
+ }
+ }
+ itk::Matrix<double,3,3> MeasurementFrame;
+ MeasurementFrame.SetIdentity();
+
+ // check ImageOrientationPatient and figure out slice direction in
+ // L-P-I (right-handed) system.
+ // In Dicom, the coordinate frame is L-P by default. Look at
+ // http://medical.nema.org/dicom/2007/07_03pu.pdf , page 301
+ itk::Matrix<double,3,3> LPSDirCos;
+
+ {
+ double dirCosArray[6];
+ // 0020,0037 -- Image Orientation (Patient)
+ allHeaders[0]->GetDirCosArray(dirCosArray);
+ double *dirCosArrayP = dirCosArray;
+ for(unsigned i = 0; i < 2; ++i)
+ {
+ for(unsigned j = 0; j < 3; ++j,++dirCosArrayP)
+ {
+ LPSDirCos[j][i] = *dirCosArrayP;
+ }
+ }
+ }
+
+ // Cross product, this gives I-axis direction
+ LPSDirCos[0][2] = (LPSDirCos[1][0]*LPSDirCos[2][1]-LPSDirCos[2][0]*LPSDirCos[1][1]);
+ LPSDirCos[1][2] = (LPSDirCos[2][0]*LPSDirCos[0][1]-LPSDirCos[0][0]*LPSDirCos[2][1]);
+ LPSDirCos[2][2] = (LPSDirCos[0][0]*LPSDirCos[1][1]-LPSDirCos[1][0]*LPSDirCos[0][1]);
+
+ std::cout << "ImageOrientationPatient (0020:0037): ";
+ std::cout << "LPS Orientation Matrix" << std::endl;
+ std::cout << LPSDirCos << std::endl;
+
+ itk::Matrix<double,3,3> SpacingMatrix;
+ SpacingMatrix.Fill(0.0);
+ SpacingMatrix[0][0]=xRes;
+ SpacingMatrix[1][1]=yRes;
+ SpacingMatrix[2][2]=sliceSpacing;
+ std::cout << "SpacingMatrix" << std::endl;
+ std::cout << SpacingMatrix << std::endl;
+
+ itk::Matrix<double,3,3> OrientationMatrix;
+ OrientationMatrix.SetIdentity();
+
+ itk::Matrix<double,3,3> NRRDSpaceDirection;
+ std::string nrrdSpaceDefinition="left-posterior-superior";;
+ NRRDSpaceDirection=LPSDirCos*OrientationMatrix*SpacingMatrix;
+
+ std::cout << "NRRDSpaceDirection" << std::endl;
+ std::cout << NRRDSpaceDirection << std::endl;
+
+ unsigned int mMosaic = 0; // number of raws in each mosaic block;
+ unsigned int nMosaic = 0; // number of columns in each mosaic block
+ unsigned int nSliceInVolume = 0;
+ unsigned int nVolume = 0;
+ bool SliceOrderIS(true);
+
+ // figure out slice order and mosaic arrangement.
+ if ( StringContains(vendor,"GE") || (StringContains(vendor,"SIEMENS") && !SliceMosaic) )
+ {
+ if(StringContains(vendor,"GE"))
+ {
+ MeasurementFrame=LPSDirCos;
+ }
+ else //SIEMENS data assumes a measurement frame that is the identity matrix.
+ {
+ MeasurementFrame.SetIdentity();
+ }
+ // has the measurement frame represented as an identity matrix.
+ double image0Origin[3];
+ allHeaders[0]->GetElementDS(0x0020, 0x0032, 3, image0Origin);
+ std::cout << "Slice 0: " << image0Origin[0] << " " << image0Origin[1] << " " << image0Origin[2] << std::endl;
+
+ // assume volume interleaving, i.e. the second dicom file stores
+ // the second slice in the same volume as the first dicom file
+ double image1Origin[3];
+ allHeaders[1]->GetElementDS(0x0020, 0x0032, 3, image1Origin);
+ std::cout << "Slice 0: " << image1Origin[0] << " " << image1Origin[1] << " " << image1Origin[2] << std::endl;
+
+ image1Origin[0] -= image0Origin[0];
+ image1Origin[1] -= image0Origin[1];
+ image1Origin[2] -= image0Origin[2];
+ double x1 = image1Origin[0]*(NRRDSpaceDirection[0][2]) +
+ image1Origin[1]*(NRRDSpaceDirection[1][2]) +
+ image1Origin[2]*(NRRDSpaceDirection[2][2]);
+ if (x1 < 0)
+ {
+ SliceOrderIS = false;
+ }
+ }
+ else if ( StringContains(vendor,"SIEMENS") && SliceMosaic )
+ {
+ MeasurementFrame.SetIdentity(); //The DICOM version of SIEMENS that uses private tags
+ // has the measurement frame represented as an identity matrix.
+ std::cout << "Siemens SliceMosaic......" << std::endl;
+
+ SliceOrderIS = false;
+
+ // for siemens mosaic image, figure out mosaic slice order from 0029|1010
+ // copy information stored in 0029,1010 into a string for parsing
+ std::string tag;
+ allHeaders[0]->GetElementOB(0x0029,0x1010, tag);
+ // parse SliceNormalVector from 0029,1010 tag
+ std::vector<double> valueArray(0);
+ int nItems = ExtractSiemensDiffusionInformation(tag, "SliceNormalVector", valueArray);
+ if (nItems != 3) // did not find enough information
+ {
+ std::cout << "Warning: Cannot find complete information on SliceNormalVector in 0029|1010" << std::endl;
+ std::cout << " Slice order may be wrong." << std::endl;
+ }
+ else if (valueArray[2] > 0)
+ {
+ SliceOrderIS = true;
+ }
+
+ // parse NumberOfImagesInMosaic from 0029,1010 tag
+ valueArray.resize(0);
+ nItems = ExtractSiemensDiffusionInformation(tag, "NumberOfImagesInMosaic", valueArray);
+ if (nItems == 0) // did not find enough information
+ {
+ std::cout << "Warning: Cannot find complete information on NumberOfImagesInMosaic in 0029|1010" << std:: endl;
+ std::cout << " Resulting image may contain empty slices." << std::endl;
+ }
+ else
+ {
+ nSliceInVolume = static_cast<int>(valueArray[0]);
+ mMosaic = static_cast<int> (ceil(sqrt(valueArray[0])));
+ nMosaic = mMosaic;
+ }
+ std::cout << "Mosaic in " << mMosaic << " X " << nMosaic
+ << " blocks (total number of blocks = " << valueArray[0] << ")." << std::endl;
+ }
+ else if (!multiSliceVolume && StringContains(vendor,"PHILIPS") && nSlice > 1)
+ // so this is not a philips multi-frame single dicom file
+ {
+ MeasurementFrame=LPSDirCos; //Philips oblique scans list the gradients with respect to the ImagePatientOrientation.
+ SliceOrderIS = true;
+
+ nSliceInVolume = numberOfSlicesPerVolume;
+ nVolume = nSlice/nSliceInVolume;
+
+ double image0Origin[3];
+ allHeaders[0]->GetElementDS(0x0020, 0x0032, 3, image0Origin);
+ std::cout << "Slice 0: " << image0Origin[0] << " " << image0Origin[1] << " " << image0Origin[2] << std::endl;
+
+ // assume volume interleaving, i.e. the second dicom file stores
+ // the second slice in the same volume as the first dicom file
+ double image1Origin[3];
+ allHeaders[nVolume]->GetElementDS(0x0020, 0x0032, 3, image1Origin);
+ std::cout << "Slice " << nVolume << ": " << image1Origin[0] << " "
+ << image1Origin[1] << " " << image1Origin[2] << std::endl;
+
+ image1Origin[0] -= image0Origin[0];
+ image1Origin[1] -= image0Origin[1];
+ image1Origin[2] -= image0Origin[2];
+ double x0 = image1Origin[0] * (NRRDSpaceDirection[0][2]) +
+ image1Origin[1] * (NRRDSpaceDirection[1][2]) +
+ image1Origin[2] * (NRRDSpaceDirection[2][2]);
+
+ // VAM - This needs more investigation -
+ // Should we default to false and change based on slice order
+ if (x0 < 0)
+ {
+ SliceOrderIS = false;
+ }
+ }
+ else if ( StringContains(vendor,"PHILIPS") && nSlice == 1)
+ {
+ // special handling for philips multi-frame dicom later.
+ }
+ else
+ {
+ std::cout << " Warning: vendor type not valid" << std::endl;
+ // treate the dicom series as an ordinary image and write a straight nrrd file.
+ WriteVolume<VolumeType>( readerOutput, outputVolumeHeaderName );
+ FreeHeaders(allHeaders);
+ return EXIT_SUCCESS;
+ }
+
+ if ( SliceOrderIS )
+ {
+ std::cout << "Slice order is IS" << std::endl;
+ }
+ else
+ {
+ std::cout << "Slice order is SI" << std::endl;
+ NRRDSpaceDirection[0][2] = -NRRDSpaceDirection[0][2];
+ NRRDSpaceDirection[1][2] = -NRRDSpaceDirection[1][2];
+ NRRDSpaceDirection[2][2] = -NRRDSpaceDirection[2][2];
+ }
+
+ std::cout << "Row: " << (NRRDSpaceDirection[0][0]) << ", "
+ << (NRRDSpaceDirection[1][0]) << ", "
+ << (NRRDSpaceDirection[2][0]) << std::endl;
+ std::cout << "Col: " << (NRRDSpaceDirection[0][1])
+ << ", " << (NRRDSpaceDirection[1][1])
+ << ", " << (NRRDSpaceDirection[2][1]) << std::endl;
+ std::cout << "Sli: " << (NRRDSpaceDirection[0][2])
+ << ", " << (NRRDSpaceDirection[1][2]) << ", "
+ << (NRRDSpaceDirection[2][2]) << std::endl;
+
+ const float orthoSliceSpacing = fabs((NRRDSpaceDirection[2][2]));
+
+ int nIgnoreVolume = 0; // Used for Philips Trace like images
+ std::vector<int> useVolume;
+
+ std::vector<float> bValues(0);
+ float maxBvalue = 0;
+ int nBaseline = 0;
+
+ // UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem is only of debug purposes.
+ std::vector< vnl_vector_fixed<double, 3> > DiffusionVectors;
+ std::vector< vnl_vector_fixed<double, 3> > UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem;
+ std::vector<int> ignorePhilipsSliceMultiFrame;
+
+ ////////////////////////////////////////////////////////////
+ // vendor dependent tags.
+ // read in gradient vectors and determin nBaseline and nMeasurement
+
+ if ( StringContains(vendor,"GE") )
+ {
+ //
+ // don't even try to convert DTI 6 Direction files
+ std::string seriesDescription;
+#if 0
+ if(allHeaders[0]->GetElementLO(0x0008,0x103e,seriesDescription,false) == EXIT_SUCCESS &&
+ StringContains(seriesDescription,"6 Directions"))
+ {
+ std::cerr << "Can't recover B-value & diffusion directions from DTI - 6 Directions scans" << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+#endif
+ nSliceInVolume = numberOfSlicesPerVolume;
+ nVolume = nSlice/nSliceInVolume;
+
+ // assume volume interleaving
+ std::cout << "Number of Slices: " << nSlice << std::endl;
+ std::cout << "Number of Volume: " << nVolume << std::endl;
+ std::cout << "Number of Slices in each volume: " << nSliceInVolume << std::endl;
+
+ for (unsigned int k = 0; k < nSlice; k += nSliceInVolume)
+ {
+ // parsing bvalue and gradient directions
+ vnl_vector_fixed<double, 3> vect3d;
+ vect3d.fill( 0 );
+ // for some weird reason this item in the GE dicom
+ // header is stored as an IS (Integer String) element.
+ ::itk::int32_t intb;
+ allHeaders[k]->GetElementISorOB(0x0043, 0x1039, intb);
+ float b = static_cast<float>(intb);
+
+ // allHeaders[k]->GetElementDS(0x0019, 0x10bb, 1, &vect3d[0]);
+ // allHeaders[k]->GetElementDS(0x0019, 0x10bc, 1, &vect3d[1]);
+ // allHeaders[k]->GetElementDS(0x0019, 0x10bd, 1, &vect3d[2]);
+ for(unsigned elementNum = 0x10bb; elementNum <= 0x10bd; ++elementNum)
+ {
+ int vecI(elementNum - 0x10bb);
+ allHeaders[k]->GetElementDSorOB(0x0019, elementNum, vect3d[vecI]);
+ }
+
+ vect3d[0] = -vect3d[0];
+ vect3d[1] = -vect3d[1];
+
+ bValues.push_back( b );
+ if (b == 0)
+ {
+ vect3d.fill( 0 );
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ DiffusionVectors.push_back(vect3d);
+ }
+ else
+ {
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ // vect3d.normalize();
+ DiffusionVectors.push_back(vect3d);
+ }
+
+ std::cout << "B-value: " << b <<
+ "; diffusion direction: " << vect3d[0] << ", " << vect3d[1] << ", " << vect3d[2] << std::endl;
+ }
+ }
+ else if ( StringContains(vendor,"PHILIPS"))
+ {
+ if(nSlice > 1 )
+ {
+ // assume volume interleaving
+ std::cout << "Number of Slices: " << nSlice << std::endl;
+ std::cout << "Number of Volumes: " << nVolume << std::endl;
+ std::cout << "Number of Slices in each volume: " << nSliceInVolume << std::endl;
+
+ std::string tmpString = "";
+ //NOTE: Philips interleaves the directions, so the all gradient directions can be
+ //determined in the first "nVolume" slices which represents the first slice from each
+ //of the gradient volumes.
+ for (unsigned int k = 0; k < nVolume; ++k)
+ {
+ std::string DiffusionDirectionality;
+ bool useSupplement49Definitions(false);
+ if(allHeaders[k]->GetElementCSorOB(0x0018,0x9075,DiffusionDirectionality,false) == EXIT_SUCCESS)
+ {
+ useSupplement49Definitions = true;
+ }
+
+ bool B0FieldFound = false;
+ double b=0.0;
+ if (useSupplement49Definitions == true )
+ {
+ B0FieldFound = allHeaders[k]->GetElementFD(0x0018,0x9087,b,false) == EXIT_SUCCESS;
+ }
+ else
+ {
+ float floatB;
+ if(allHeaders[k]->GetElementFLorOB(0x2001,0x1003,floatB,false) == EXIT_SUCCESS)
+ {
+ B0FieldFound = true;
+ }
+ if(B0FieldFound)
+ {
+ b = static_cast<double>(floatB);
+ }
+ std::string tag;
+ allHeaders[k]->GetElementCSorOB(0x2001, 0x1004, tag, false );
+ if(StringContains(tag,"I") && b != 0)
+ {
+ DiffusionDirectionality="ISOTROPIC";
+ }
+ }
+
+ vnl_vector_fixed<double, 3> vect3d;
+ vect3d.fill( 0 );
+ //std::cout << "HACK: " << "DiffusionDirectionality=" << DiffusionDirectionality << ", k= " << k << std::endl;
+ //std::cout << "HACK: " << "B0FieldFound=" << B0FieldFound << ", b=" << b << ", DiffusionDirectionality=" << DiffusionDirectionality << std::endl;
+
+ if ( StringContains(DiffusionDirectionality,"ISOTROPIC") )
+ { //Deal with images that are to be ignored
+ //std::cout << " SKIPPING ISOTROPIC Diffusion. " << std::endl;
+ //std::cout << "HACK: IGNORE IMAGEFILE: " << k << " of " << filenames.size() << " " << filenames[k] << std::endl;
+ // Ignore the Trace like image
+ ++nIgnoreVolume;
+ useVolume.push_back(0);
+ continue;
+ }
+ else if (( !B0FieldFound || b == 0 ) || StringContains(DiffusionDirectionality, "NONE") )
+ { //Deal with b0 images
+ bValues.push_back(b);
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ DiffusionVectors.push_back(vect3d);
+ useVolume.push_back(1);
+ continue;
+ }
+ else if ( StringContains(DiffusionDirectionality,"DIRECTIONAL") || ( DiffusionDirectionality == "" ))
+ { //Deal with gradient direction images
+ bValues.push_back(b);
+ //std::cout << "HACK: GRADIENT IMAGEFILE: " << k << " of " << filenames.size() << " " << filenames[k] << std::endl;
+ useVolume.push_back(1);
+ if (useSupplement49Definitions == true )
+ {
+ double doubleArray[3];
+ //Use alternate method to get value out of a sequence header (Some Phillips Data).
+ if(allHeaders[k]->GetElementFD(0x0018, 0x9089, 3, doubleArray, false) != EXIT_SUCCESS)
+ {
+ //std::cout << "Looking for 0018|9089 in sequence 0018,9076" << std::endl;
+ // gdcm::SeqEntry *
+ // DiffusionSeqEntry=allHeaders[k]->GetSeqEntry(0x0018,0x9076);
+ itk::DCMTKSequence DiffusionSeqEntry;
+ allHeaders[k]->GetElementSQ(0x0018,0x9076,DiffusionSeqEntry);
+ // const unsigned int
+ // n=DiffusionSeqEntry->GetNumberOfSQItems();
+ unsigned int n = DiffusionSeqEntry.card();
+ if( n == 0 )
+ {
+ std::cout << "ERROR: Sequence entry 0018|9076 has no items." << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+ DiffusionSeqEntry.GetElementFD(0x0018,0x9089, 3, doubleArray);
+ }
+ vect3d[0] = doubleArray[0];
+ vect3d[1] = doubleArray[1];
+ vect3d[2] = doubleArray[2];
+ std::cout << "===== gradient orientations:" << k << " "
+ << inputFileNames[k] << " (0018,9089) " << " " << vect3d << std::endl;
+ }
+ else
+ {
+ float tmp[3];
+ /*const bool b0exist =*/
+ allHeaders[k]->GetElementFLorOB( 0x2005, 0x10b0, tmp[0] );
+ allHeaders[k]->GetElementFLorOB( 0x2005, 0x10b1, tmp[1] );
+ allHeaders[k]->GetElementFLorOB( 0x2005, 0x10b2, tmp[2] );
+ vect3d[0] = static_cast<double>(tmp[0]);
+ vect3d[1] = static_cast<double>(tmp[1]);
+ vect3d[2] = static_cast<double>(tmp[2]);
+ }
+
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ // vect3d.normalize();
+ DiffusionVectors.push_back(vect3d);
+ }
+ else // Have no idea why we'd be here so error out
+ {
+ std::cout << "ERROR: DiffusionDirectionality was "
+ << DiffusionDirectionality << " Don't know what to do with that..." << std::endl;
+ FreeHeaders(allHeaders);
+ return EXIT_FAILURE;
+ }
+
+ std::cout << "B-value: " << b <<
+ "; diffusion direction: " << vect3d[0] << ", " << vect3d[1] << ", " << vect3d[2] << std::endl;
+ }
+ }
+ else
+ {
+ // multi-frame file, everything is inside
+ std::map<std::vector<double>, double> gradientDirectionAndBValue;
+ ignorePhilipsSliceMultiFrame.clear();
+
+ sliceLocations.clear();
+ bValues.clear();
+ DiffusionVectors.clear();
+ useVolume.clear();
+
+ itk::DCMTKSequence perFrameFunctionalGroup;
+ itk::DCMTKSequence innerSeq;
+ double dwbValue;
+
+ allHeaders[0]->GetElementSQ(0x5200,0x9230,perFrameFunctionalGroup);
+ int nItems = perFrameFunctionalGroup.card();
+
+ // have to determine if volume slices are interleaved
+ std::string origins[2];
+
+ for(unsigned long i = 0;
+ i < static_cast<unsigned long>(perFrameFunctionalGroup.card()); ++i)
+ {
+ itk::DCMTKItem curItem;
+ perFrameFunctionalGroup.GetElementItem(i,curItem);
+
+ // index slice locations with string origin
+ itk::DCMTKSequence originSeq;
+ curItem.GetElementSQ(0x0020, 0x9113,originSeq);
+ std::string originString;
+ originSeq.GetElementDS(0x0020,0x0032,originString);
+ ++sliceLocations[originString];
+ // save origin of first 2 slices to compare and see if the
+ // volume is interleaved.
+ if(i < 2)
+ {
+ origins[i] = originString;
+ }
+
+ itk::DCMTKSequence mrDiffusionSeq;
+ curItem.GetElementSQ(0x0018,0x9117,mrDiffusionSeq);
+
+ std::string dirValue;
+ mrDiffusionSeq.GetElementCSorOB(0x0018,0x9075,dirValue);
+
+ if ( StringContains(dirValue,"ISO") )
+ {
+ useVolume.push_back(0);
+ ignorePhilipsSliceMultiFrame.push_back( i );
+ }
+ else if ( StringContains(dirValue,"NONE") )
+ {
+ useVolume.push_back(1);
+ std::vector<double> v(3);
+ v[0] = 0; v[1] = 0; v[2] = 0;
+ unsigned int nOld = gradientDirectionAndBValue.size();
+ gradientDirectionAndBValue[v] = 0;
+ unsigned int nNew = gradientDirectionAndBValue.size();
+
+ if (nOld != nNew)
+ {
+ vnl_vector_fixed<double, 3> vect3d;
+ vect3d.fill( 0 );
+ DiffusionVectors.push_back( vect3d );
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ bValues.push_back( 0 );
+
+ }
+ }
+ else
+ {
+ useVolume.push_back(1);
+ if(mrDiffusionSeq.GetElementDSorOB(0x0018,0x9087,dwbValue,false) != EXIT_SUCCESS)
+ {
+ mrDiffusionSeq.GetElementFD(0x0018,0x9087,dwbValue);
+ }
+ itk::DCMTKSequence volSeq;
+ mrDiffusionSeq.GetElementSQ(0x0018, 0x9076,volSeq);
+ double dwgVal[3];
+ if(volSeq.GetElementDSorOB<double>(0x0018,0x9089,3,dwgVal,false) != EXIT_SUCCESS)
+ {
+ volSeq.GetElementFD(0x0018,0x9089,3,dwgVal);
+ }
+ std::vector<double> v(3);
+ v[0] = dwgVal[0];
+ v[1] = dwgVal[1];
+ v[2] = dwgVal[2];
+ unsigned int nOld = gradientDirectionAndBValue.size();
+ gradientDirectionAndBValue[v] = dwbValue;
+ unsigned int nNew = gradientDirectionAndBValue.size();
+
+ if (nOld != nNew)
+ {
+ vnl_vector_fixed<double, 3> vect3d;
+ vect3d[0] = v[0]; vect3d[1] = v[1]; vect3d[2] = v[2];
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ // vect3d.normalize();
+ DiffusionVectors.push_back( vect3d );
+
+ bValues.push_back( dwbValue);
+ }
+ }
+ }
+
+ numberOfSlicesPerVolume=sliceLocations.size();
+
+ // de-interleave slices if the origins of the first 2 slices
+ // are the same.
+ if(origins[0] == origins[1])
+ {
+ // interleaved image
+ DeInterleaveVolume(readerOutput,numberOfSlicesPerVolume,perFrameFunctionalGroup.card());
+ }
+
+
+ std::cout << "LPS Matrix: " << std::endl << LPSDirCos << std::endl;
+ std::cout << "Volume Origin: " << std::endl << ImageOrigin[0] << ","
+ << ImageOrigin[1] << "," << ImageOrigin[2] << "," << std::endl;
+ std::cout << "Number of slices per volume: " << numberOfSlicesPerVolume << std::endl;
+ std::cout << "Slice matrix size: " << nRows << " X " << nCols << std::endl;
+ std::cout << "Image resolution: " << xRes << ", " << yRes << ", " << sliceSpacing << std::endl;
+
+ NRRDSpaceDirection=LPSDirCos*OrientationMatrix*SpacingMatrix;
+
+ MeasurementFrame=LPSDirCos;
+
+ nSliceInVolume = sliceLocations.size();
+ nVolume = nItems/nSliceInVolume;
+ nIgnoreVolume = ignorePhilipsSliceMultiFrame.size()/nSliceInVolume;
+
+ for( unsigned int k2 = 0; k2 < bValues.size(); ++k2 )
+ {
+ std::cout << k2 << ": direction: " << DiffusionVectors[k2][0]
+ << ", " << DiffusionVectors[k2][1] << ", " << DiffusionVectors[k2][2]
+ << ", b-value: " << bValues[k2] << std::endl;
+ }
+
+ }
+ }
+ else if ( StringContains(vendor, "SIEMENS") )
+ {
+
+ int nStride = 1;
+
+ if ( !SliceMosaic )
+ {
+ std::cout << orthoSliceSpacing << std::endl;
+ nSliceInVolume = numberOfSlicesPerVolume;
+ nVolume = nSlice/nSliceInVolume;
+ std::cout << "Number of Slices: " << nSlice << std::endl;
+ std::cout << "Number of Volume: " << nVolume << std::endl;
+ std::cout << "Number of Slices in each volume: " << nSliceInVolume << std::endl;
+ nStride = nSliceInVolume;
+ }
+ else
+ {
+ std::cout << "Data in Siemens Mosaic Format" << std::endl;
+ nVolume = nSlice;
+ std::cout << "Number of Volume: " << nVolume << std::endl;
+ std::cout << "Number of Slices in each volume: " << nSliceInVolume << std::endl;
+ nStride = 1;
+ }
+
+ // JTM - Determine bvalues from all gradients
+ vnl_vector_fixed<double, 3> vect3d;
+
+ for (unsigned int k = 0; k < nSlice; k += nStride )
+ {
+ //
+ // in Siemens, this entry is a 'CSA Header' which is blob
+ // of mixed text & binary data. Pretty annoying but there you
+ // have it.
+ std::string diffusionInfoString;;
+ allHeaders[k]->GetElementOB( 0x0029, 0x1010, diffusionInfoString );
+
+ // parse B_value from 0029,1010 tag
+ std::vector<double> valueArray(0);
+
+ int nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", valueArray);
+ if (nItems != 1)
+ {
+ //
+ // B_Value is missing -- the punt position is to count this
+ // volume as having a B_value & Gradient Direction of zero
+ std::cout << "Warning: Cannot find complete information on B_value in 0029|1010" << std::endl;
+ bValues.push_back( 0.0 );
+ vect3d.fill( 0.0 );
+ UnmodifiedDiffusionVectorsInDicomLPSCoordinateSystem.push_back(vect3d);
+ DiffusionVectors.push_back(vect3d);
+ continue;
+ }
+
+ // we got a 'valid' B-value
+ // If we're trusting the gradient directions in the header,
+ // then all we need to do here is save the bValue.
+ if(!useBMatrixGradientDirections)
+ {
+ valueArray.resize(0);
+ ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", valueArray);
+
+ bValues.push_back( valueArray[0] );
+ }
+ else
+ {
+ // JTM - Patch from UNC: fill the nhdr header with the gradient directions and
+ // bvalues computed out of the BMatrix
+ valueArray.resize(0);
+ nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_matrix", valueArray);
+ vnl_matrix_fixed<double, 3, 3> bMatrix;