From b4c0426aa9e3a0bed7b2c5af46fc6c809e50f02b Mon Sep 17 00:00:00 2001 From: chaircrusher Date: Wed, 4 Sep 2013 09:15:19 -0500 Subject: [PATCH 1/2] BUG: fix generation of vectors from CSA data --- DWIConvert/SiemensDWIConverter.h | 169 +++++++++++++++++++++++++++++++ 1 file changed, 169 insertions(+) diff --git a/DWIConvert/SiemensDWIConverter.h b/DWIConvert/SiemensDWIConverter.h index f96fe063b6..b4c611bce6 100644 --- a/DWIConvert/SiemensDWIConverter.h +++ b/DWIConvert/SiemensDWIConverter.h @@ -20,6 +20,175 @@ class SiemensDWIConverter : public DWIConverter } virtual ~SiemensDWIConverter() {} +<<<<<<< HEAD +======= + template + T CSAExtractFromString(const char *ptr) + { + T rval = *(reinterpret_cast(ptr)); + itk::ByteSwapper::SwapFromSystemToLittleEndian(&rval); + return rval; + } + + class CSAItem : public std::vector + { + public: + typedef std::vector SuperClass; + + itk::uint32_t vm; + std::string vr; + + CSAItem(unsigned int length) : SuperClass(length), + vm(0) + { + } + CSAItem() : SuperClass() + { + } + CSAItem( const CSAItem & other ) : SuperClass(other.size()) + { + *this = other; + } + CSAItem & operator=(const CSAItem &other) + { + this->resize(0); + for(CSAItem::const_iterator it = other.begin(); + it != other.end(); + ++it) + { + this->push_back(*it); + } + this->vm = other.vm; + this->vr = other.vr; + return *this; + } + + template + std::vector AsVector() const + { + std::vector rval; + for(unsigned i = 0; i < this->size(); ++i) + { + if(! (*this)[i].empty()) + { + T val = 0; + std::stringstream convert((*this)[i]); + convert >> val; + rval.push_back(val); + } + } + return rval; + } + void DebugPrint() const + { + std::cerr << " VM = " << this->vm << " VR = " << this->vr << std::endl + << " "; + bool firstTime(false); + for(CSAItem::const_iterator it = this->begin(); + it != this->end(); ++it) + { + if(firstTime) + { + firstTime = false; + } + else + { + std::cerr << " "; + } + std::cerr << *it; + } + std::cerr << std::endl; + } + }; + + class CSAHeader : public std::map + { + public: + void DebugPrint() const + { + for(CSAHeader::const_iterator it = this->begin(); + it != this->end(); ++it) + { + std::cerr << it->first << std::endl; + it->second.DebugPrint(); + } + } + }; + + void DecodeCSAHeader(CSAHeader &header, const std::string &infoString) + { + // + // the reference used to write this code is here: + // http://nipy.sourceforge.net/nibabel/dicom/siemens_csa.html + const char *info = infoString.c_str(); + + const bool isCSA2 = info[0] == 'S' && info[1] == 'V' + && info[2] == '1' && info[3] == '0'; + unsigned int offset; + + if(isCSA2) + { + offset = 8; // past SV10 + unused 4 bytes + } + else + { + offset = 0; + } + const itk::uint32_t numberOfTags = + this->CSAExtractFromString(info+offset); + offset += sizeof(itk::uint32_t); // skip numberOfTags; + offset += sizeof(itk::uint32_t); // skip unused2 + + for(unsigned i = 0; i < numberOfTags; ++i) + { + // tag name is 64 bytes null terminated. + std::string tagName = info + offset; + offset += 64; // skip tag name + itk::uint32_t vm = this->CSAExtractFromString(info+offset); + offset += sizeof(itk::uint32_t); + + CSAItem current(vm); + current.vm = vm; + + // vr = 3 bytes of string + 1 for pad + char vr[4]; + for(unsigned j = 0; j < 3; ++j) + { + vr[j] = info[offset + j]; + } + vr[3] = '\0'; + current.vr = vr; + offset += 4; // after VR + offset += 4; // skip syngodt + + const itk::int32_t nItems = + this->CSAExtractFromString(info + offset); + offset += 4; + offset += 4; // skip xx + + for(int j = 0; j < nItems; ++j) + { + // 4 items in XX, first being item length + const itk::int32_t itemLength = + this->CSAExtractFromString(info + offset); + offset += 16; + std::string valueString; + valueString = info + offset; + offset += itemLength; + while((offset % 4) != 0) + { + ++offset; + } + if(j < static_cast(vm)) + { + current[j] = valueString; + } + } + header[tagName] = current; + } + } + +>>>>>>> BUG: fix generation of vectors from CSA data /** Siemens datasets are either in the * normal sequential volume arrangement or * in mosaic datasets, where each slice contains From 19962ae21b065be2fe9d1c252cdb1b5c6bcb259f Mon Sep 17 00:00:00 2001 From: chaircrusher Date: Wed, 4 Sep 2013 11:21:11 -0500 Subject: [PATCH 2/2] BUG: core dump from referencing wrong variable --- DWIConvert/SiemensDWIConverter.h | 92 ++++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 29 deletions(-) diff --git a/DWIConvert/SiemensDWIConverter.h b/DWIConvert/SiemensDWIConverter.h index b4c611bce6..9a947051e7 100644 --- a/DWIConvert/SiemensDWIConverter.h +++ b/DWIConvert/SiemensDWIConverter.h @@ -20,8 +20,6 @@ class SiemensDWIConverter : public DWIConverter } virtual ~SiemensDWIConverter() {} -<<<<<<< HEAD -======= template T CSAExtractFromString(const char *ptr) { @@ -188,7 +186,6 @@ class SiemensDWIConverter : public DWIConverter } } ->>>>>>> BUG: fix generation of vectors from CSA data /** Siemens datasets are either in the * normal sequential volume arrangement or * in mosaic datasets, where each slice contains @@ -268,16 +265,28 @@ class SiemensDWIConverter : public DWIConverter std::string diffusionInfoString;; this->m_Headers[k]->GetElementOB( 0x0029, 0x1010, diffusionInfoString ); + CSAHeader csaHeader; + this->DecodeCSAHeader(csaHeader,diffusionInfoString); + // csaHeader.DebugPrint(); + // parse B_value from 0029,1010 tag std::vector valueArray(0); - // 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. + CSAHeader::const_iterator csaIt; if( !this->m_UseBMatrixGradientDirections ) { - valueArray.resize(0); - int nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", valueArray); - + int nItems(0); + if((csaIt = csaHeader.find("B_value")) == csaHeader.end()) + { + nItems = 0; + } + else + { + // 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. + valueArray = csaIt->second.AsVector(); + nItems = valueArray.size(); + } vnl_vector_fixed vect3d; if( nItems != 1 ) @@ -291,6 +300,15 @@ class SiemensDWIConverter : public DWIConverter continue; } this->m_BValues.push_back( valueArray[0] ); + if((csaIt = csaHeader.find("DiffusionGradientDirection")) == csaHeader.end()) + { + nItems = 0; + } + else + { + valueArray = csaIt->second.AsVector(); + nItems = valueArray.size(); + } // parse DiffusionGradientDirection from 0029,1010 tag valueArray.resize(0); @@ -340,29 +358,48 @@ class SiemensDWIConverter : public DWIConverter } else { + int nItems; + // UNC comments: We get the value of the b-value tag in the header. + // We won't use it as is, but just to locate the B0 images. + // This check must be added, otherwise the bmatrix of the B0 is not + // read properly (it's not an actual field in the DICOM header of the B0). + std::vector bval_tmp(0); + bool b0_image = false; + // UNC comments: Get the bvalue + if((csaIt = csaHeader.find("B_value")) == csaHeader.end()) + { + nItems = 0; + } + else + { + // 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. + bval_tmp = csaIt->second.AsVector(); + } + + if( bval_tmp[0] == 0 ) + { + b0_image = true; + } + // JTM - Patch from UNC: fill the nhdr header with the gradient directions and // bvalues computed out of the BMatrix - valueArray.resize(0); - int nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_matrix", valueArray); + if((csaIt = csaHeader.find("B_matrix")) == csaHeader.end()) + { + nItems = 0; + } + else + { + valueArray = csaIt->second.AsVector(); + nItems = valueArray.size(); + } vnl_matrix_fixed bMatrix; - if( nItems == 6 ) + if( nItems == 6 && !b0_image) { std::cout << "=============================================" << std::endl; std::cout << "BMatrix calculations..." << std::endl; - // UNC comments: We get the value of the b-value tag in the header. - // We won't use it as is, but just to locate the B0 images. - // This check must be added, otherwise the bmatrix of the B0 is not - // read properly (it's not an actual field in the DICOM header of the B0). - std::vector bval_tmp(0); - bool b0_image = false; - - // UNC comments: Get the bvalue - nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", bval_tmp); - if( bval_tmp[0] == 0 ) - { - b0_image = true; - } // UNC comments: The principal eigenvector of the bmatrix is to be extracted as // it's the gradient direction and trace of the matrix is the b-value @@ -420,10 +457,7 @@ class SiemensDWIConverter : public DWIConverter } else { - // silently returning zero gradient vectors is a problem, - // but it is also necessary for some files. - valueArray.resize(0); - ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", valueArray); + this->m_BValues.push_back( bval_tmp[0] ); vnl_vector_fixed vect3d; this->m_BValues.push_back( valueArray[0] ); vect3d[0] = 0;