#include <iostream>
int main( int argc, char *argv[] )
{
if(argc < 3)
{
std::cerr << "Usage: " << argv[0] << " NrrdFileName(.nhdr) threshold(on B0)"
<< " FAImageFileName RelativeAnisotropyFileName " <<
"[ExtractGradientAndReferenceImage from the NRRD file and "
<< "write them as images]" << std::endl;
std::cerr << "\tExample args: xt_dwi.nhdr 80 FA.mhd 1" << std::endl;
return EXIT_FAILURE;
}
const unsigned int Dimension = 3;
unsigned int numberOfImages = 0;
unsigned int numberOfGradientImages = 0;
bool readb0 = false;
double b0 = 0;
typedef unsigned short PixelType;
ImageType::Pointer img;
reader->SetFileName(argv[1]);
try
{
reader->Update();
img = reader->GetOutput();
}
{
std::cout << ex << std::endl;
return EXIT_FAILURE;
}
PixelType, PixelType, double > TensorReconstructionImageFilterType;
std::vector<std::string> imgMetaKeys = imgMetaDictionary.
GetKeys();
std::vector<std::string>::const_iterator itKey = imgMetaKeys.begin();
std::string metaString;
TensorReconstructionImageFilterType::GradientDirectionType vect3d;
TensorReconstructionImageFilterType::GradientDirectionContainerType::Pointer
DiffusionVectors =
TensorReconstructionImageFilterType::GradientDirectionContainerType::New();
for (; itKey != imgMetaKeys.end(); ++itKey)
{
double x,y,z;
itk::ExposeMetaData<std::string> (imgMetaDictionary, *itKey, metaString);
if (itKey->find("DWMRI_gradient") != std::string::npos)
{
std::cout << *itKey << " ---> " << metaString << std::endl;
sscanf(metaString.c_str(), "%lf %lf %lf\n", &x, &y, &z);
vect3d[0] = x; vect3d[1] = y; vect3d[2] = z;
DiffusionVectors->InsertElement( numberOfImages, vect3d );
++numberOfImages;
if (vect3d[0] == 0.0 &&
vect3d[1] == 0.0 &&
vect3d[2] == 0.0)
{
continue;
}
++numberOfGradientImages;
}
else if (itKey->find("DWMRI_b-value") != std::string::npos)
{
std::cout << *itKey << " ---> " << metaString << std::endl;
readb0 = true;
b0 = atof(metaString.c_str());
}
}
std::cout << "Number of gradient images: "
<< numberOfGradientImages
<< " and Number of reference images: "
<< numberOfImages - numberOfGradientImages
<< std::endl;
if(!readb0)
{
std::cerr << "BValue not specified in header file" << std::endl;
return EXIT_FAILURE;
}
typedef ReferenceImageType GradientImageType;
std::vector< GradientImageType::Pointer > imageContainer;
DWIIteratorType dwiit( img, img->GetBufferedRegion() );
for( unsigned int i = 0; i<numberOfImages; i++ )
{
GradientImageType::Pointer image = GradientImageType::New();
image->CopyInformation( img );
image->SetBufferedRegion( img->GetBufferedRegion() );
image->SetRequestedRegion( img->GetRequestedRegion() );
image->Allocate();
IteratorType it( image, image->GetBufferedRegion() );
dwiit.GoToBegin();
it.GoToBegin();
while (!it.IsAtEnd())
{
it.Set(dwiit.Get()[i]);
++it;
++dwiit;
}
imageContainer.push_back( image );
}
if( (argc > 4) && (atoi(argv[6])) )
{
unsigned int referenceImageIndex = 0;
for( unsigned int i = 0; i<numberOfImages; i++ )
{
GradientWriterType::Pointer gradientWriter = GradientWriterType::New();
gradientWriter->SetInput( imageContainer[i] );
char filename[100];
if (DiffusionVectors->ElementAt(i).two_norm() <= 0.0)
{
std::string fn("ReferenceImage%d.mhd");
sprintf(filename, fn.c_str(), referenceImageIndex );
++referenceImageIndex;
}
else
{
std::string fn("Gradient%d.mhd");
sprintf(filename, fn.c_str(), i );
}
gradientWriter->SetFileName( filename );
gradientWriter->Update();
}
}
TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter =
TensorReconstructionImageFilterType::New();
tensorReconstructionFilter->SetGradientImage( DiffusionVectors, reader->GetOutput() );
tensorReconstructionFilter->SetNumberOfThreads( 1 );
tensorReconstructionFilter->SetBValue(b0);
tensorReconstructionFilter->SetThreshold( static_cast<
TensorReconstructionImageFilterType::ReferencePixelType >(
atof(argv[2])));
tensorReconstructionFilter->Update();
TensorReconstructionImageFilterType::OutputImageType > TensorWriterType;
TensorWriterType::Pointer tensorWriter = TensorWriterType::New();
tensorWriter->SetFileName( argv[3] );
tensorWriter->SetInput( tensorReconstructionFilter->GetOutput() );
tensorWriter->Update();
typedef TensorReconstructionImageFilterType::TensorPixelType TensorPixelType;
typedef TensorPixelType::RealValueType RealValueType;
TensorReconstructionImageFilterType::OutputImageType, FAImageType >
FAFilterType;
FAFilterType::Pointer fractionalAnisotropyFilter = FAFilterType::New();
fractionalAnisotropyFilter->SetInput( tensorReconstructionFilter->GetOutput() );
FAWriterType::Pointer faWriter = FAWriterType::New();
faWriter->SetInput( fractionalAnisotropyFilter->GetOutput() );
faWriter->SetFileName(argv[4]);
faWriter->Update();
typedef TensorReconstructionImageFilterType::TensorPixelType TensorPixelType;
typedef TensorPixelType::RealValueType RealValueType;
TensorReconstructionImageFilterType::OutputImageType, RAImageType >
RAFilterType;
RAFilterType::Pointer relativeAnisotropyFilter = RAFilterType::New();
relativeAnisotropyFilter->SetInput( tensorReconstructionFilter->GetOutput() );
RAWriterType::Pointer raWriter = RAWriterType::New();
raWriter->SetInput( relativeAnisotropyFilter->GetOutput() );
raWriter->SetFileName(argv[5]);
raWriter->Update();
return EXIT_SUCCESS;
}