Skip to content

Commit

Permalink
COMP: Undo experimental changes in DWIConvert
Browse files Browse the repository at this point in the history
  • Loading branch information
chaircrusher committed Sep 3, 2013
1 parent bdcd136 commit 41a5b34
Showing 1 changed file with 15 additions and 231 deletions.
246 changes: 15 additions & 231 deletions DWIConvert/SiemensDWIConverter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,6 @@
#define __SiemensDWIConverter_h
#include "DWIConverter.h"
#include "StringContains.h"
#include "itkIntTypes.h"
#include "itkByteSwapper.h"
#include <vector>
#include <sstream>

/** specific converter for Siemens scanners*/
class SiemensDWIConverter : public DWIConverter
Expand All @@ -24,172 +20,6 @@ class SiemensDWIConverter : public DWIConverter
}
virtual ~SiemensDWIConverter() {}

template <typename T>
T CSAExtractFromString(const char *ptr)
{
T rval = *(reinterpret_cast<const T *>(ptr));
itk::ByteSwapper<T>::SwapFromSystemToLittleEndian(&rval);
return rval;
}

class CSAItem : public std::vector<std::string >
{
public:
typedef std::vector<std::string> 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 <typename T>
std::vector<T> AsVector() const
{
std::vector<T> rval(this->size());
for(unsigned i = 0; i < this->size(); ++i)
{
T val = 0;
if(! (*this)[i].empty())
{
std::stringstream convert((*this)[i]);
convert >> val;
}
rval[i] = 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<std::string,CSAItem>
{
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<itk::uint32_t>(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<itk::uint32_t>(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<itk::int32_t>(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<itk::int32_t>(info + offset);
offset += 16;
std::string valueString;
valueString = info + offset;
offset += itemLength;
while((offset % 4) != 0)
{
++offset;
}
if(j < static_cast<int>(vm))
{
current[j] = valueString;
}
}
header[tagName] = current;
}
}

/** Siemens datasets are either in the
* normal sequential volume arrangement or
* in mosaic datasets, where each slice contains
Expand Down Expand Up @@ -269,28 +99,16 @@ 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<double> 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.
if( !this->m_UseBMatrixGradientDirections )
{
int nItems(0);
CSAHeader::const_iterator it = csaHeader.find("B_value");
if(it == 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 = it->second.AsVector<double>();
nItems = valueArray.size();
}
valueArray.resize(0);
int nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", valueArray);

vnl_vector_fixed<double, 3> vect3d;

if( nItems != 1 )
Expand All @@ -304,16 +122,10 @@ class SiemensDWIConverter : public DWIConverter
continue;
}
this->m_BValues.push_back( valueArray[0] );
if(csaHeader.find("DiffusionGradientDirection") == csaHeader.end())
{
nItems = 0;
}
else
{
valueArray = csaHeader["DiffusionGradientDirection"].AsVector<double>();
nItems = valueArray.size();
}

// parse DiffusionGradientDirection from 0029,1010 tag
valueArray.resize(0);
nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "DiffusionGradientDirection", valueArray);
if( nItems != 3 ) // did not find enough information
{
std::cout << "Warning: Cannot find complete information on DiffusionGradientDirection in 0029|1010"
Expand Down Expand Up @@ -361,16 +173,8 @@ class SiemensDWIConverter : public DWIConverter
{
// JTM - Patch from UNC: fill the nhdr header with the gradient directions and
// bvalues computed out of the BMatrix
int nItems;
if(csaHeader.find("B_matrix") == csaHeader.end())
{
nItems = 0;
}
else
{
valueArray = csaHeader["B_matrix"].AsVector<double>();
nItems = valueArray.size();
}
valueArray.resize(0);
int nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_matrix", valueArray);
vnl_matrix_fixed<double, 3, 3> bMatrix;

if( nItems == 6 )
Expand All @@ -385,18 +189,7 @@ class SiemensDWIConverter : public DWIConverter
bool b0_image = false;

// UNC comments: Get the bvalue
if(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 = csaHeader["B_Value"].AsVector<double>();
}

nItems = ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", bval_tmp);
if( bval_tmp[0] == 0 )
{
b0_image = true;
Expand Down Expand Up @@ -460,19 +253,10 @@ class SiemensDWIConverter : public DWIConverter
{
// silently returning zero gradient vectors is a problem,
// but it is also necessary for some files.
if(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 = csaHeader["B_Value"].AsVector<double>();
}
this->m_BValues.push_back( valueArray[0] );
valueArray.resize(0);
ExtractSiemensDiffusionInformation(diffusionInfoString, "B_value", valueArray);
vnl_vector_fixed<double, 3> vect3d;
this->m_BValues.push_back( valueArray[0] );
vect3d[0] = 0;
vect3d[1] = 0;
vect3d[2] = 0;
Expand Down

0 comments on commit 41a5b34

Please sign in to comment.