mitk: mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 Copyright (c) German Cancer Research Center, Division of Medical and
00003 Biological Informatics. All rights reserved.
00004 See MITKCopyright.txt or https://www.mitk.org/copyright.html for details.
00005 
00006 This software is distributed WITHOUT ANY WARRANTY; without even
00007 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00008 PURPOSE.  See the above copyright notices for more information.
00009 
00010 =========================================================================*/
00011 
00012 
00013 #include "mitkSiemensMosaicDicomDiffusionImageHeaderReader.h"
00014 
00015 #include "gdcmGlobal.h"
00016 //#include "gdcmVersion.h"
00017 
00018 #if  GDCM_MAJOR_VERSION >= 2
00019 #define DGDCM2
00020 #endif
00021 
00022 #ifndef DGDCM2
00023 
00024 #include "gdcm.h"
00025 
00026 // relevant Siemens private tags
00027 const gdcm::DictEntry SiemensMosiacParameters( 0x0051, 0x100b, "IS", "1", "Mosiac Matrix Size" );
00028 const gdcm::DictEntry SiemensDictNMosiac( 0x0019, 0x100a, "US", "1", "Number of Images In Mosaic" );
00029 const gdcm::DictEntry SiemensDictBValue( 0x0019, 0x100c, "IS", "1", "B Value of diffusion weighting" );
00030 const gdcm::DictEntry SiemensDictDiffusionDirection( 0x0019, 0x100e, "FD", "3", "Diffusion Gradient Direction" );
00031 const gdcm::DictEntry SiemensDictDiffusionMatrix( 0x0019, 0x1027, "FD", "6", "Diffusion Matrix" );
00032 const gdcm::DictEntry SiemensDictShadowInfo( 0x0029, 0x1010, "OB", "1", "Siemens DWI Info" );
00033 
00034 #else
00035 
00036 #include "gdcmFile.h"
00037 #include "gdcmImageReader.h"
00038 #include "gdcmDict.h"
00039 #include "gdcmDicts.h"
00040 #include "gdcmDictEntry.h"
00041 #include "gdcmDictEntry.h"
00042 #include "gdcmDict.h"
00043 #include "gdcmFile.h"
00044 #include "gdcmSerieHelper.h"
00045 
00046 #endif
00047 
00048 
00049 mitk::SiemensMosaicDicomDiffusionImageHeaderReader::SiemensMosaicDicomDiffusionImageHeaderReader()
00050 {
00051 }
00052 
00053 mitk::SiemensMosaicDicomDiffusionImageHeaderReader::~SiemensMosaicDicomDiffusionImageHeaderReader()
00054 {
00055 }
00056 
00057 int mitk::SiemensMosaicDicomDiffusionImageHeaderReader::ExtractSiemensDiffusionInformation( std::string tagString, std::string nameString, std::vector<double>& valueArray )
00058 {
00059   std::string::size_type atPosition = tagString.find( nameString );
00060   if ( atPosition == std::string::npos)
00061   {
00062     return 0;
00063   }
00064   else
00065   {
00066     std::string infoAsString = tagString.substr( atPosition, tagString.size()-atPosition+1 );
00067     const char * infoAsCharPtr = infoAsString.c_str();
00068 
00069     int vm = *(infoAsCharPtr+64);
00070     std::string vr = infoAsString.substr( 68, 4 );
00071     //int syngodt = *(infoAsCharPtr+72);
00072     //int nItems = *(infoAsCharPtr+76);
00073     //int localDummy = *(infoAsCharPtr+80);
00074 
00075     int offset = 84;
00076     for (int k = 0; k < vm; k++)
00077     {
00078       int itemLength = *(infoAsCharPtr+offset+4);
00079       int strideSize = static_cast<int> (ceil(static_cast<double>(itemLength)/4) * 4);
00080       std::string valueString = infoAsString.substr( offset+16, itemLength );
00081       valueArray.push_back( atof(valueString.c_str()) );
00082       offset += 16+strideSize;
00083     }
00084     return vm;
00085   }
00086 }
00087 
00088 // do the work
00089 void mitk::SiemensMosaicDicomDiffusionImageHeaderReader::Update()
00090 {
00091 
00092   // check if there are filenames
00093   if(m_DicomFilenames.size())
00094   {
00095     // adapted from namic-sandbox
00096     // DicomToNrrdConverter.cxx
00097 
00098     VolumeReaderType::DictionaryArrayRawPointer inputDict
00099         = m_VolumeReader->GetMetaDataDictionaryArray();
00100 
00101 #ifndef DGDCM2
00102     if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensMosiacParameters.GetKey()) == 0)
00103       gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensMosiacParameters);
00104     if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictNMosiac.GetKey()) == 0)
00105       gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictNMosiac);
00106     if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictBValue.GetKey()) == 0)
00107       gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictBValue);
00108     if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictDiffusionDirection.GetKey()) == 0)
00109       gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictDiffusionDirection);
00110     if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictDiffusionMatrix.GetKey()) == 0)
00111       gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictDiffusionMatrix);
00112     if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictShadowInfo.GetKey()) == 0)
00113       gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictShadowInfo);
00114 #endif
00115 
00116     ReadPublicTags();
00117 
00118     int mMosaic = 0;   // number of raws in each mosaic block;
00119     int nMosaic = 0;   // number of columns in each mosaic block
00120 
00121     std::cout << "Siemens SliceMosaic......" << std::endl;
00122 
00123     m_SliceOrderIS = false;
00124 
00125     // for siemens mosaic image, figure out mosaic slice order from 0029|1010
00126     std::string tag;
00127     tag.clear();
00128 
00129 #ifndef DGDCM2
00130 
00131     gdcm::File *header0 = new gdcm::File;
00132     gdcm::BinEntry* binEntry;
00133 
00134     header0->SetMaxSizeLoadEntry(65536);
00135     header0->SetFileName( m_DicomFilenames[0] );
00136     header0->SetLoadMode( gdcm::LD_ALL );
00137     header0->Load();
00138 
00139     // copy information stored in 0029,1010 into a string for parsing
00140     gdcm::DocEntry* docEntry = header0->GetFirstEntry();
00141     while(docEntry)
00142     {
00143       if ( docEntry->GetKey() == "0029|1010"  )
00144       {
00145         binEntry = dynamic_cast<gdcm::BinEntry*> ( docEntry );
00146         int binLength = binEntry->GetFullLength();
00147         tag.resize( binLength );
00148         uint8_t * tagString = binEntry->GetBinArea();
00149 
00150         for (int n = 0; n < binLength; n++)
00151         {
00152           tag[n] = *(tagString+n);
00153         }
00154         break;
00155       }
00156       docEntry = header0->GetNextEntry();
00157     }
00158 
00159 #else
00160 
00161     gdcm::ImageReader reader;
00162     reader.SetFileName( m_DicomFilenames[0].c_str() );
00163     if( !reader.Read() )
00164     {
00165       itkExceptionMacro(<< "Cannot read requested file");
00166     }
00167     const gdcm::File &f = reader.GetFile();
00168     const gdcm::DataSet &ds = f.GetDataSet();
00169 
00170     //    gdcm::DataSet ds = header0->GetDataSet();
00171     gdcm::DataSet::ConstIterator it = ds.Begin();
00172 
00173     // Copy of the header->content
00174     // copy information stored in 0029,1010 into a string for parsing
00175     for(; it != ds.End(); ++it)
00176     {
00177       const gdcm::DataElement &ref = *it;
00178       if (ref.GetTag() == gdcm::Tag(0x0029,0x1010)) {
00179         tag = std::string(ref.GetByteValue()->GetPointer(),ref.GetByteValue()->GetLength());
00180       }
00181     }
00182 
00183 #endif
00184 
00185     // parse SliceNormalVector from 0029,1010 tag
00186     std::vector<double> valueArray(0);
00187     int nItems = ExtractSiemensDiffusionInformation(tag, "SliceNormalVector", valueArray);
00188     if (nItems != 3)  // did not find enough information
00189     {
00190       std::cout << "Warning: Cannot find complete information on SliceNormalVector in 0029|1010\n";
00191       std::cout << "         Slice order may be wrong.\n";
00192     }
00193     else if (valueArray[2] > 0)
00194     {
00195       m_SliceOrderIS = true;
00196     }
00197 
00198     // parse NumberOfImagesInMosaic from 0029,1010 tag
00199     valueArray.resize(0);
00200     nItems = ExtractSiemensDiffusionInformation(tag, "NumberOfImagesInMosaic", valueArray);
00201     if (nItems == 0)  // did not find enough information
00202     {
00203       std::cout << "Warning: Cannot find complete information on NumberOfImagesInMosaic in 0029|1010\n";
00204       std::cout << "         Resulting image may contain empty slices.\n";
00205     }
00206     else 
00207     {
00208       this->m_Output->nSliceInVolume = static_cast<int>(valueArray[0]);
00209       mMosaic = static_cast<int> (ceil(sqrt(valueArray[0])));
00210       nMosaic = mMosaic;
00211     }
00212     std::cout << "Mosaic in " << mMosaic << " X " << nMosaic << " blocks (total number of blocks = " << valueArray[0] << ").\n"; 
00213 
00214 
00215     ReadPublicTags2();
00216 
00217     int nStride = 1;
00218 
00219     std::cout << "Data in Siemens Mosaic Format\n";
00220     //nVolume = nSlice;
00221     //std::cout << "Number of Volume: " << nVolume << std::endl;
00222     std::cout << "Number of Slices in each volume: " << this->m_Output->nSliceInVolume << std::endl;
00223     nStride = 1;
00224 
00225     for (int k = 0; k < m_nSlice; k += nStride )
00226     {
00227 
00228 #ifndef DGDCM2
00229 
00230       gdcm::File *header0 = new gdcm::File;
00231       gdcm::BinEntry* binEntry;
00232 
00233       header0->SetMaxSizeLoadEntry(65536);
00234       header0->SetFileName( m_DicomFilenames[0] );
00235       header0->SetLoadMode( gdcm::LD_ALL );
00236       header0->Load();
00237 
00238       // copy information stored in 0029,1010 into a string for parsing
00239       gdcm::DocEntry* docEntry = header0->GetFirstEntry();
00240       while(docEntry)
00241       {
00242         if ( docEntry->GetKey() == "0029|1010"  )
00243         {
00244           binEntry = dynamic_cast<gdcm::BinEntry*> ( docEntry );
00245           int binLength = binEntry->GetFullLength();
00246           tag.resize( binLength );
00247           uint8_t * tagString = binEntry->GetBinArea();
00248 
00249           for (int n = 0; n < binLength; n++)
00250           {
00251             tag[n] = *(tagString+n);
00252           }
00253           break;
00254         }
00255         docEntry = header0->GetNextEntry();
00256       }
00257 
00258 #else
00259 
00260       gdcm::ImageReader reader;
00261       reader.SetFileName( m_DicomFilenames[0].c_str() );
00262       if( !reader.Read() )
00263       {
00264         itkExceptionMacro(<< "Cannot read requested file");
00265       }
00266       const gdcm::File &f = reader.GetFile();
00267       const gdcm::DataSet &ds = f.GetDataSet();
00268 
00269       //    gdcm::DataSet ds = header0->GetDataSet();
00270       gdcm::DataSet::ConstIterator it = ds.Begin();
00271 
00272       // Copy of the header->content
00273       // copy information stored in 0029,1010 into a string for parsing
00274       for(; it != ds.End(); ++it)
00275       {
00276         const gdcm::DataElement &ref = *it;
00277         if (ref.GetTag() == gdcm::Tag(0x0029,0x1010)) {
00278           tag = std::string(ref.GetByteValue()->GetPointer(),ref.GetByteValue()->GetLength());
00279         }
00280       }
00281 
00282 #endif
00283 
00284       // parse B_value from 0029,1010 tag
00285       std::vector<double> valueArray(0);
00286       vnl_vector_fixed<double, 3> vect3d;
00287       int nItems = ExtractSiemensDiffusionInformation(tag, "B_value", valueArray);
00288       if (nItems != 1 || valueArray[0] == 0)  // did not find enough information
00289       {
00290         MITK_INFO << "Reading diffusion info from 0019|100c and 0019|100e tags";
00291         tag.clear();
00292         bool success = itk::ExposeMetaData<std::string> ( *(*inputDict)[0], "0019|100c", tag );
00293         if(success)
00294         {
00295           this->m_Output->bValue = atof( tag.c_str() );
00296           tag.clear();
00297           success = itk::ExposeMetaData<std::string> ( *(*inputDict)[k], "0019|100e",  tag);
00298           if(success)
00299           {
00300             memcpy( &vect3d[0], tag.c_str()+0, 8 );
00301             memcpy( &vect3d[1], tag.c_str()+8, 8 );
00302             memcpy( &vect3d[2], tag.c_str()+16, 8 );
00303             vect3d.normalize();
00304             this->m_Output->DiffusionVector = vect3d;
00305             TransformGradients();
00306             MITK_INFO << "BV: " << this->m_Output->bValue;
00307             MITK_INFO << " GD: " << this->m_Output->DiffusionVector << std::endl;
00308             continue;
00309           }
00310         }
00311       }
00312       else
00313       {
00314         MITK_INFO << "Reading diffusion info from 0029|1010 tags";
00315         this->m_Output->bValue = valueArray[0];
00316 
00317         // parse DiffusionGradientDirection from 0029,1010 tag
00318         valueArray.resize(0);
00319         nItems = ExtractSiemensDiffusionInformation(tag, "DiffusionGradientDirection", valueArray);
00320         if (nItems == 3)
00321         {
00322           vect3d[0] = valueArray[0];
00323           vect3d[1] = valueArray[1];
00324           vect3d[2] = valueArray[2];
00325           vect3d.normalize();
00326           this->m_Output->DiffusionVector = vect3d;
00327           TransformGradients();
00328           MITK_INFO << "BV: " << this->m_Output->bValue;
00329           MITK_INFO << " GD: " << this->m_Output->DiffusionVector;
00330           continue;
00331         }
00332       }
00333 
00334       MITK_ERROR << "No diffusion info found, assuming BASELINE" << std::endl;
00335       this->m_Output->bValue = 0.0;
00336       vect3d.fill( 0.0 );
00337       this->m_Output->DiffusionVector = vect3d;
00338 
00339     }
00340 
00341   }
00342 }