diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp index 0e5753d1e6..507f3c183a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp @@ -1,533 +1,533 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkNrrdDiffusionImageReader_cpp #define __mitkNrrdDiffusionImageReader_cpp #include "mitkNrrdDiffusionImageReader.h" #include "itkImageFileReader.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkNiftiImageIO.h" #include #include #include "itksys/SystemTools.hxx" namespace mitk { template void NrrdDiffusionImageReader ::GenerateData() { // Since everything is completely read in GenerateOutputInformation() it is stored // in a cache variable. A timestamp is associated. // If the timestamp of the cache variable is newer than the MTime, we only need to // assign the cache variable to the DataObject. // Otherwise, the tree must be read again from the file and OuputInformation must // be updated! if ( ( ! m_OutputCache ) || ( this->GetMTime( ) > m_CacheTime.GetMTime( ) ) ) { this->GenerateOutputInformation(); itkWarningMacro("Cache regenerated!"); } if (!m_OutputCache) { itkWarningMacro("cache is empty!"); } static_cast(this->GetOutput()) ->SetVectorImage(m_OutputCache->GetVectorImage()); static_cast(this->GetOutput()) ->SetB_Value(m_OutputCache->GetB_Value()); static_cast(this->GetOutput()) ->SetDirections(m_OutputCache->GetDirections()); static_cast(this->GetOutput()) ->SetMeasurementFrame(m_OutputCache->GetMeasurementFrame()); static_cast(this->GetOutput()) ->InitializeFromVectorImage(); } template void NrrdDiffusionImageReader::GenerateOutputInformation() { typename OutputType::Pointer outputForCache = OutputType::New(); if ( m_FileName == "") { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename to be read is empty!"); } else { try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } MITK_INFO << "NrrdDiffusionImageReader: reading image information"; typename ImageType::Pointer img; std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".hdwi" || ext == ".dwi") { typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(this->m_FileName); itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); reader->SetImageIO(io); reader->Update(); img = reader->GetOutput(); int vecsize = img->GetVectorLength(); std::cout << vecsize << std::endl; } else if(ext == ".fsl" || ext == ".fslgz") { // create temporary file with correct ending for nifti-io - std::string fname3 = m_FileName; + std::string fname3 = "temp_dwi"; fname3 += ext == ".fsl" ? ".nii" : ".nii.gz"; itksys::SystemTools::CopyAFile(m_FileName.c_str(), fname3.c_str()); // create reader and read file typedef itk::Image ImageType4D; itk::NiftiImageIO::Pointer io2 = itk::NiftiImageIO::New(); typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(fname3); reader->SetImageIO(io2); reader->Update(); typename ImageType4D::Pointer img4 = reader->GetOutput(); // delete temporary file itksys::SystemTools::RemoveFile(fname3.c_str()); // convert 4D file to vector image img = ImageType::New(); typename ImageType::SpacingType spacing; typename ImageType4D::SpacingType spacing4 = img4->GetSpacing(); for(int i=0; i<3; i++) spacing[i] = spacing4[i]; img->SetSpacing( spacing ); // Set the image spacing typename ImageType::PointType origin; typename ImageType4D::PointType origin4 = img4->GetOrigin(); for(int i=0; i<3; i++) origin[i] = origin4[i]; img->SetOrigin( origin ); // Set the image origin typename ImageType::DirectionType direction; typename ImageType4D::DirectionType direction4 = img4->GetDirection(); for(int i=0; i<3; i++) for(int j=0; j<3; j++) direction[i][j] = direction4[i][j]; img->SetDirection( direction ); // Set the image direction typename ImageType::RegionType region; typename ImageType4D::RegionType region4 = img4->GetLargestPossibleRegion(); typename ImageType::RegionType::SizeType size; typename ImageType4D::RegionType::SizeType size4 = region4.GetSize(); for(int i=0; i<3; i++) size[i] = size4[i]; typename ImageType::RegionType::IndexType index; typename ImageType4D::RegionType::IndexType index4 = region4.GetIndex(); for(int i=0; i<3; i++) index[i] = index4[i]; region.SetSize(size); region.SetIndex(index); img->SetRegions( region ); img->SetVectorLength(size4[3]); img->Allocate(); itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); typedef typename ImageType::PixelType VecPixType; for (it = it.Begin(); !it.IsAtEnd(); ++it) { VecPixType vec = it.Get(); typename ImageType::IndexType currentIndex = it.GetIndex(); for(int i=0; i<3; i++) index4[i] = currentIndex[i]; for(unsigned int ind=0; indGetPixel(index4); } it.Set(vec); } } m_DiffusionVectors = GradientDirectionContainerType::New(); m_OriginalDiffusionVectors = GradientDirectionContainerType::New(); if (ext == ".hdwi" || ext == ".dwi") { itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); std::vector::const_iterator itKey = imgMetaKeys.begin(); std::string metaString; GradientDirectionType vect3d; int numberOfImages = 0; int numberOfGradientImages = 0; bool readb0 = false; bool readFrame = false; double xx, xy, xz, yx, yy, yz, zx, zy, zz; for (; itKey != imgMetaKeys.end(); itKey ++) { double x,y,z; itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); if (itKey->find("DWMRI_gradient") != std::string::npos) { sscanf(metaString.c_str(), "%lf %lf %lf\n", &x, &y, &z); vect3d[0] = x; vect3d[1] = y; vect3d[2] = z; m_DiffusionVectors->InsertElement( numberOfImages, vect3d ); ++numberOfImages; // If the direction is 0.0, this is a reference image 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) { readb0 = true; m_B_Value = atof(metaString.c_str()); } else if (itKey->find("measurement frame") != std::string::npos) { sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz); readFrame = true; if (xx>10e-10 || xy>10e-10 || xz>10e-10 || yx>10e-10 || yy>10e-10 || yz>10e-10 || zx>10e-10 || zy>10e-10 || zz>10e-10 ) { m_MeasurementFrame(0,0) = xx; m_MeasurementFrame(0,1) = xy; m_MeasurementFrame(0,2) = xz; m_MeasurementFrame(1,0) = yx; m_MeasurementFrame(1,1) = yy; m_MeasurementFrame(1,2) = yz; m_MeasurementFrame(2,0) = zx; m_MeasurementFrame(2,1) = zy; m_MeasurementFrame(2,2) = zz; } else { m_MeasurementFrame(0,0) = 1; m_MeasurementFrame(0,1) = 0; m_MeasurementFrame(0,2) = 0; m_MeasurementFrame(1,0) = 0; m_MeasurementFrame(1,1) = 1; m_MeasurementFrame(1,2) = 0; m_MeasurementFrame(2,0) = 0; m_MeasurementFrame(2,1) = 0; m_MeasurementFrame(2,2) = 1; } } } if(!readb0) { MITK_INFO << "BValue not specified in header file"; } } else if(ext == ".fsl" || ext == ".fslgz") { std::string line; std::vector bvec_entries; std::string fname = m_FileName; fname += ".bvecs"; std::ifstream myfile (fname.c_str()); if (myfile.is_open()) { while ( myfile.good() ) { getline (myfile,line); char* pch = strtok (const_cast(line.c_str())," "); while (pch != NULL) { bvec_entries.push_back(atof(pch)); pch = strtok (NULL, " "); } } myfile.close(); } else { MITK_INFO << "Unable to open bvecs file"; } std::vector bval_entries; std::string fname2 = m_FileName; fname2 += ".bvals"; std::ifstream myfile2 (fname2.c_str()); if (myfile2.is_open()) { while ( myfile2.good() ) { getline (myfile2,line); char* pch = strtok (const_cast(line.c_str())," "); while (pch != NULL) { bval_entries.push_back(atof(pch)); pch = strtok (NULL, " "); } } myfile2.close(); } else { MITK_INFO << "Unable to open bvals file"; } m_B_Value = -1; unsigned int numb = bval_entries.size(); for(unsigned int i=0; i vec; vec[0] = bvec_entries.at(i); vec[1] = bvec_entries.at(i+numb); vec[2] = bvec_entries.at(i+2*numb); // Adjust the vector length to encode gradient strength float factor = b_val/m_B_Value; if(vec.magnitude() > 0) { vec[0] = sqrt(factor)*vec[0]; vec[1] = sqrt(factor)*vec[1]; vec[2] = sqrt(factor)*vec[2]; } m_DiffusionVectors->InsertElement(i,vec); } for(int i=0; i<3; i++) for(int j=0; j<3; j++) m_MeasurementFrame[i][j] = i==j ? 1 : 0; } outputForCache->SetVectorImage(img); outputForCache->SetB_Value(m_B_Value); outputForCache->SetDirections(m_DiffusionVectors); outputForCache->SetMeasurementFrame(m_MeasurementFrame); // Since we have already read the tree, we can store it in a cache variable // so that it can be assigned to the DataObject in GenerateData(); m_OutputCache = outputForCache; m_CacheTime.Modified(); try { setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } } catch(std::exception& e) { MITK_INFO << "Std::Exception while reading file!!"; MITK_INFO << e.what(); throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); } catch(...) { MITK_INFO << "Exception while reading file!!"; throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested vessel tree file!"); } } } template const char* NrrdDiffusionImageReader ::GetFileName() const { return m_FileName.c_str(); } template void NrrdDiffusionImageReader ::SetFileName(const char* aFileName) { m_FileName = aFileName; } template const char* NrrdDiffusionImageReader ::GetFilePrefix() const { return m_FilePrefix.c_str(); } template void NrrdDiffusionImageReader ::SetFilePrefix(const char* aFilePrefix) { m_FilePrefix = aFilePrefix; } template const char* NrrdDiffusionImageReader ::GetFilePattern() const { return m_FilePattern.c_str(); } template void NrrdDiffusionImageReader ::SetFilePattern(const char* aFilePattern) { m_FilePattern = aFilePattern; } template bool NrrdDiffusionImageReader ::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/) { // First check the extension if( filename == "" ) { return false; } std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".hdwi" || ext == ".dwi") { itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetImageIO(io); reader->SetFileName(filename); try { reader->Update(); } catch(itk::ExceptionObject e) { MITK_INFO << e.GetDescription(); } typename ImageType::Pointer img = reader->GetOutput(); itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); std::vector::const_iterator itKey = imgMetaKeys.begin(); std::string metaString; for (; itKey != imgMetaKeys.end(); itKey ++) { itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); if (itKey->find("modality") != std::string::npos) { if (metaString.find("DWMRI") != std::string::npos) { return true; } } } } if (ext == ".fsl" || ext == ".fslgz") { // itk::NiftiImageIO::Pointer io2 = itk::NiftiImageIO::New(); // typedef itk::ImageFileReader FileReaderType; // typename FileReaderType::Pointer reader = FileReaderType::New(); // reader->SetImageIO(io2); // reader->SetFileName(filename); // try // { // reader->Update(); // } // catch(itk::ExceptionObject e) // { // MITK_INFO << e.GetDescription(); // } std::string fname = filename; fname += ".bvecs"; std::string fname2 = filename; fname2 += ".bvals"; if( itksys::SystemTools::FileExists(fname.c_str()) && itksys::SystemTools::FileExists(fname2.c_str()) ) { return true; } else { MITK_INFO << ".bvals and .bvals files do not exist properly"; } } return false; } } //namespace MITK #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp index 6a7a366e84..7ab9f2ec1b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp @@ -1,326 +1,459 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkNrrdTensorImageReader.h" #include "itkImageFileReader.h" #include "itkImageRegionIterator.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" - +#include #include "mitkITKImageImport.h" #include "mitkImageDataItem.h" namespace mitk { void NrrdTensorImageReader ::GenerateData() { if ( m_FileName == "") { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!"); } else { try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } - MITK_INFO << "Loading tensor image ..."; - - typedef itk::VectorImage ImageType; - itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); - typedef itk::ImageFileReader FileReaderType; - FileReaderType::Pointer reader = FileReaderType::New(); - reader->SetImageIO(io); - reader->SetFileName(this->m_FileName); - reader->Update(); - ImageType::Pointer img = reader->GetOutput(); - - typedef itk::Image,3> VecImgType; - VecImgType::Pointer vecImg = VecImgType::New(); - vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing - vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin - vecImg->SetDirection( img->GetDirection() ); // Set the image direction - vecImg->SetRegions( img->GetLargestPossibleRegion()); - vecImg->Allocate(); - - itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); - ot = ot.Begin(); - - itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); - it = it.Begin(); - - typedef ImageType::PixelType VarPixType; - typedef VecImgType::PixelType FixPixType; - int numComponents = img->GetNumberOfComponentsPerPixel(); - - itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); - std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); - std::vector::const_iterator itKey = imgMetaKeys.begin(); - std::string metaString; - - bool readFrame = false; - double xx, xy, xz, yx, yy, yz, zx, zy, zz; - MeasurementFrameType measFrame; - measFrame.SetIdentity(); - MeasurementFrameType measFrameTransp; - measFrameTransp.SetIdentity(); - - for (; itKey != imgMetaKeys.end(); itKey ++) + try { - itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); - if (itKey->find("measurement frame") != std::string::npos) - { - sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz); + MITK_INFO << "Trying to load dti as nifti ..."; + std::string fname3 = "temp_dti.nii"; + itksys::SystemTools::CopyAFile(m_FileName.c_str(), fname3.c_str()); - if (xx>10e-10 || xy>10e-10 || xz>10e-10 || - yx>10e-10 || yy>10e-10 || yz>10e-10 || - zx>10e-10 || zy>10e-10 || zz>10e-10 ) - { - readFrame = true; - - measFrame(0,0) = xx; - measFrame(0,1) = xy; - measFrame(0,2) = xz; - measFrame(1,0) = yx; - measFrame(1,1) = yy; - measFrame(1,2) = yz; - measFrame(2,0) = zx; - measFrame(2,1) = zy; - measFrame(2,2) = zz; - - measFrameTransp = measFrame.GetTranspose(); - } - } - } - - if (numComponents==6) - { - while (!it.IsAtEnd()) - { - // T'=RTR' - VarPixType vec = it.Get(); - FixPixType fixVec(vec.GetDataPointer()); + typedef itk::Image ImageType; + itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New(); + typedef itk::ImageFileReader FileReaderType; + FileReaderType::Pointer reader = FileReaderType::New(); + reader->SetImageIO(io); + reader->SetFileName(fname3); + reader->Update(); + itksys::SystemTools::RemoveFile(fname3.c_str()); - if(readFrame) - { - itk::DiffusionTensor3D tensor; - tensor.SetElement(0, vec.GetElement(0)); - tensor.SetElement(1, vec.GetElement(1)); - tensor.SetElement(2, vec.GetElement(2)); - tensor.SetElement(3, vec.GetElement(3)); - tensor.SetElement(4, vec.GetElement(4)); - tensor.SetElement(5, vec.GetElement(5)); + ImageType::Pointer img = reader->GetOutput(); + itk::Size<4> size = img->GetLargestPossibleRegion().GetSize(); - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); - fixVec = tensor; - } + itk::ImageRegion< 3 > region; + region.SetSize(0,size[0]); + region.SetSize(1,size[1]); + region.SetSize(2,size[2]); + itk::Vector spacing; + spacing[0] = img->GetSpacing()[0]; + spacing[1] = img->GetSpacing()[1]; + spacing[2] = img->GetSpacing()[2]; + itk::Point origin; + origin[0] = img->GetOrigin()[0]; + origin[1] = img->GetOrigin()[1]; + origin[2] = img->GetOrigin()[2]; + itk::Matrix direction; + direction[0][0] = img->GetDirection()[0][0]; + direction[1][0] = img->GetDirection()[1][0]; + direction[2][0] = img->GetDirection()[2][0]; + direction[0][1] = img->GetDirection()[0][1]; + direction[1][1] = img->GetDirection()[1][1]; + direction[2][1] = img->GetDirection()[2][1]; + direction[0][2] = img->GetDirection()[0][2]; + direction[1][2] = img->GetDirection()[1][2]; + direction[2][2] = img->GetDirection()[2][2]; + + typedef itk::Image,3> VecImgType; + typedef VecImgType::PixelType FixPixType; + VecImgType::Pointer vecImg = VecImgType::New(); + vecImg->SetSpacing( spacing ); + vecImg->SetOrigin( origin ); + vecImg->SetDirection( direction ); + vecImg->SetRegions( region ); + vecImg->Allocate(); + + itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); + ot = ot.Begin(); - ot.Set(fixVec); - ++ot; - ++it; - } - } - else if(numComponents==9) - { - while (!it.IsAtEnd()) + while (!ot.IsAtEnd()) { - VarPixType vec = it.Get(); itk::DiffusionTensor3D tensor; - tensor.SetElement(0, vec.GetElement(0)); - tensor.SetElement(1, vec.GetElement(1)); - tensor.SetElement(2, vec.GetElement(2)); - tensor.SetElement(3, vec.GetElement(4)); - tensor.SetElement(4, vec.GetElement(5)); - tensor.SetElement(5, vec.GetElement(8)); - - if(readFrame) + ImageType::IndexType idx; + idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2]; + + if (size[3]==6) + { + for (int te=0; teGetPixel(idx)); + } + } + else if (size[3]==9) { - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + idx[3] = 0; + tensor.SetElement(0, img->GetPixel(idx)); + idx[3] = 1; + tensor.SetElement(1, img->GetPixel(idx)); + idx[3] = 2; + tensor.SetElement(2, img->GetPixel(idx)); + idx[3] = 4; + tensor.SetElement(3, img->GetPixel(idx)); + idx[3] = 5; + tensor.SetElement(4, img->GetPixel(idx)); + idx[3] = 8; + tensor.SetElement(5, img->GetPixel(idx)); } + else + throw itk::ImageFileReaderException(__FILE__, __LINE__, "Unknown number of komponents for DTI file. Should be 6 or 9!"); FixPixType fixVec(tensor); ot.Set(fixVec); ++ot; - ++it; } + this->GetOutput()->InitializeByItk(vecImg.GetPointer()); + this->GetOutput()->SetVolume(vecImg->GetBufferPointer()); } - else if (numComponents==1) + catch(...) { - typedef itk::Image ImageType; + MITK_INFO << "Trying to load dti as nrrd ..."; + + typedef itk::VectorImage ImageType; itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetImageIO(io); reader->SetFileName(this->m_FileName); reader->Update(); ImageType::Pointer img = reader->GetOutput(); - itk::Size<4> size = img->GetLargestPossibleRegion().GetSize(); + typedef itk::Image,3> VecImgType; + VecImgType::Pointer vecImg = VecImgType::New(); + vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing + vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin + vecImg->SetDirection( img->GetDirection() ); // Set the image direction + vecImg->SetRegions( img->GetLargestPossibleRegion()); + vecImg->Allocate(); + + itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); + ot = ot.Begin(); + + itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); + it = it.Begin(); + + typedef ImageType::PixelType VarPixType; + typedef VecImgType::PixelType FixPixType; + int numComponents = img->GetNumberOfComponentsPerPixel(); + + itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); + std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); + std::vector::const_iterator itKey = imgMetaKeys.begin(); + std::string metaString; + + bool readFrame = false; + double xx, xy, xz, yx, yy, yz, zx, zy, zz; + MeasurementFrameType measFrame; + measFrame.SetIdentity(); + MeasurementFrameType measFrameTransp; + measFrameTransp.SetIdentity(); + + for (; itKey != imgMetaKeys.end(); itKey ++) + { + itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); + if (itKey->find("measurement frame") != std::string::npos) + { + sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz); + + if (xx>10e-10 || xy>10e-10 || xz>10e-10 || + yx>10e-10 || yy>10e-10 || yz>10e-10 || + zx>10e-10 || zy>10e-10 || zz>10e-10 ) + { + readFrame = true; + + measFrame(0,0) = xx; + measFrame(0,1) = xy; + measFrame(0,2) = xz; + measFrame(1,0) = yx; + measFrame(1,1) = yy; + measFrame(1,2) = yz; + measFrame(2,0) = zx; + measFrame(2,1) = zy; + measFrame(2,2) = zz; + + measFrameTransp = measFrame.GetTranspose(); + } + } + } - while (!ot.IsAtEnd()) + if (numComponents==6) { - itk::DiffusionTensor3D tensor; - ImageType::IndexType idx; - idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2]; - for (int te=0; teGetPixel(idx)); + // T'=RTR' + VarPixType vec = it.Get(); + FixPixType fixVec(vec.GetDataPointer()); + + if(readFrame) + { + itk::DiffusionTensor3D tensor; + tensor.SetElement(0, vec.GetElement(0)); + tensor.SetElement(1, vec.GetElement(1)); + tensor.SetElement(2, vec.GetElement(2)); + tensor.SetElement(3, vec.GetElement(3)); + tensor.SetElement(4, vec.GetElement(4)); + tensor.SetElement(5, vec.GetElement(5)); + + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + fixVec = tensor; + } + + ot.Set(fixVec); + ++ot; + ++it; } - if(readFrame) + } + else if(numComponents==9) + { + while (!it.IsAtEnd()) { - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + VarPixType vec = it.Get(); + itk::DiffusionTensor3D tensor; + tensor.SetElement(0, vec.GetElement(0)); + tensor.SetElement(1, vec.GetElement(1)); + tensor.SetElement(2, vec.GetElement(2)); + tensor.SetElement(3, vec.GetElement(4)); + tensor.SetElement(4, vec.GetElement(5)); + tensor.SetElement(5, vec.GetElement(8)); + + if(readFrame) + { + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + } + + FixPixType fixVec(tensor); + ot.Set(fixVec); + ++ot; + ++it; } - FixPixType fixVec(tensor); - ot.Set(fixVec); - ++ot; } - } - else - { - throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!"); - } + else if (numComponents==1) + { + typedef itk::Image ImageType; + itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); + typedef itk::ImageFileReader FileReaderType; + FileReaderType::Pointer reader = FileReaderType::New(); + reader->SetImageIO(io); + reader->SetFileName(this->m_FileName); + reader->Update(); + ImageType::Pointer img = reader->GetOutput(); + + itk::Size<4> size = img->GetLargestPossibleRegion().GetSize(); - this->GetOutput()->InitializeByItk(vecImg.GetPointer()); - this->GetOutput()->SetVolume(vecImg->GetBufferPointer()); + MITK_INFO << size; + + while (!ot.IsAtEnd()) + { + itk::DiffusionTensor3D tensor; + ImageType::IndexType idx; + idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2]; + + if (size[3]==6) + { + for (int te=0; teGetPixel(idx)); + } + + // idx[3] = 0; + // tensor.SetElement(0, img->GetPixel(idx)); + // idx[3] = 1; + // tensor.SetElement(1, img->GetPixel(idx)); + // idx[3] = 3; + // tensor.SetElement(2, img->GetPixel(idx)); + // idx[3] = 2; + // tensor.SetElement(3, img->GetPixel(idx)); + // idx[3] = 4; + // tensor.SetElement(4, img->GetPixel(idx)); + // idx[3] = 5; + // tensor.SetElement(5, img->GetPixel(idx)); + } + else if (size[3]==9) + { + idx[3] = 0; + tensor.SetElement(0, img->GetPixel(idx)); + idx[3] = 1; + tensor.SetElement(1, img->GetPixel(idx)); + idx[3] = 2; + tensor.SetElement(2, img->GetPixel(idx)); + idx[3] = 4; + tensor.SetElement(3, img->GetPixel(idx)); + idx[3] = 5; + tensor.SetElement(4, img->GetPixel(idx)); + idx[3] = 8; + tensor.SetElement(5, img->GetPixel(idx)); + } + else + throw itk::ImageFileReaderException(__FILE__, __LINE__, "Unknown number of komponents for DTI file. Should be 6 or 9!"); + + if(readFrame) + { + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + } + FixPixType fixVec(tensor); + ot.Set(fixVec); + ++ot; + } + } + else + { + throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!"); + } + + this->GetOutput()->InitializeByItk(vecImg.GetPointer()); + this->GetOutput()->SetVolume(vecImg->GetBufferPointer()); + + } try { setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } } catch(std::exception& e) { throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); } catch(...) { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested DTI file!"); } } } void NrrdTensorImageReader::GenerateOutputInformation() { } const char* NrrdTensorImageReader ::GetFileName() const { return m_FileName.c_str(); } void NrrdTensorImageReader ::SetFileName(const char* aFileName) { m_FileName = aFileName; } const char* NrrdTensorImageReader ::GetFilePrefix() const { return m_FilePrefix.c_str(); } void NrrdTensorImageReader ::SetFilePrefix(const char* aFilePrefix) { m_FilePrefix = aFilePrefix; } const char* NrrdTensorImageReader ::GetFilePattern() const { return m_FilePattern.c_str(); } void NrrdTensorImageReader ::SetFilePattern(const char* aFilePattern) { m_FilePattern = aFilePattern; } bool NrrdTensorImageReader ::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/) { // First check the extension if( filename == "" ) { return false; } std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".hdti" || ext == ".dti") { return true; } return false; } itk::DiffusionTensor3D NrrdTensorImageReader ::ConvertMatrixTypeToFixedArrayType(const itk::DiffusionTensor3D::Superclass::MatrixType & matrix) { /* | 0 1 2 | * | X 3 4 | * | X X 5 | */ itk::DiffusionTensor3D arr; arr.SetElement(0,matrix(0,0)); arr.SetElement(1,matrix(0,1)); arr.SetElement(2,matrix(0,2)); arr.SetElement(3,matrix(1,3)); arr.SetElement(4,matrix(1,4)); arr.SetElement(5,matrix(2,5)); return arr; } } //namespace MITK diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp b/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp index 0fe0498d0b..45cadfb53b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp @@ -1,114 +1,111 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include #include #include #include #include #include #include #include #include #include #include "itkVectorContainer.h" #include "vnl/vnl_vector_fixed.h" #include "itkVectorImage.h" #include "mitkBaseDataIOFactory.h" #include "mitkDiffusionImage.h" #include "mitkBaseData.h" #include "mitkDiffusionCoreObjectFactory.h" #include "mitkCoreObjectFactory.h" #include "mitkCoreExtObjectFactory.h" #include "mitkImageWriter.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include +#include using namespace mitk; /** * Convert files from one ending to the other */ int QballReconstruction(int argc, char* argv[]) { - if ( argc!=5 + if ( argc!=6 ) { std::cout << std::endl; std::cout << "Perform QBall reconstruction on an dwi file" << std::endl; std::cout << std::endl; - std::cout << "usage: " << argv[0] << " " << std::endl; + std::cout << "usage: " << argv[0] << " " << argv[1] << " " << std::endl; std::cout << std::endl; return EXIT_FAILURE; } try { RegisterDiffusionCoreObjectFactory(); MITK_INFO << "Loading image ..."; const std::string s1="", s2=""; - std::vector infile = BaseDataIO::LoadBaseDataFromFile( argv[1], s1, s2, false ); + std::vector infile = BaseDataIO::LoadBaseDataFromFile( argv[2], s1, s2, false ); DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; dwi->AverageRedundantGradients(0.001); FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); filter->SetBValue(dwi->GetB_Value()); - filter->SetThreshold( boost::lexical_cast(argv[3]) ); - filter->SetLambda(boost::lexical_cast(argv[2])); + filter->SetThreshold( boost::lexical_cast(argv[4]) ); + filter->SetLambda(boost::lexical_cast(argv[3])); filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); - std::string outfilename = argv[4]; + std::string outfilename = argv[5]; MITK_INFO << "writing image " << outfilename; - mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); - for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) - { - if ( (*it)->CanWriteBaseDataType(image.GetPointer()) ) { - (*it)->SetFileName( outfilename.c_str() ); - (*it)->DoWrite( image.GetPointer() ); - } - } + mitk::NrrdQBallImageWriter::Pointer writer = mitk::NrrdQBallImageWriter::New(); + writer->SetInput(image.GetPointer()); + writer->SetFileName(outfilename.c_str()); + writer->Update(); } catch ( itk::ExceptionObject &err) { MITK_INFO << "Exception: " << err; } catch ( std::exception err) { MITK_INFO << "Exception: " << err.what(); } catch ( ... ) { MITK_INFO << "Exception!"; } return EXIT_SUCCESS; } RegisterDiffusionCoreMiniApp(QballReconstruction); diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/TensorReconstruction.cpp b/Modules/DiffusionImaging/DiffusionCore/MiniApps/TensorReconstruction.cpp index fad1f84050..1f9eea61a3 100644 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/TensorReconstruction.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/TensorReconstruction.cpp @@ -1,168 +1,89 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -#include - -#include "itkVectorContainer.h" -#include "vnl/vnl_vector_fixed.h" -#include "itkVectorImage.h" - #include "mitkBaseDataIOFactory.h" #include "mitkDiffusionImage.h" #include "mitkBaseData.h" #include "mitkDiffusionCoreObjectFactory.h" -#include "mitkCoreObjectFactory.h" -#include "mitkCoreExtObjectFactory.h" -#include "mitkImageWriter.h" - -#include "itkDiffusionTensor3DReconstructionImageFilter.h" -#include "mitkTensorImage.h" +#include +#include +#include +#include +using namespace mitk; /** * Convert files from one ending to the other */ int TensorReconstruction(int argc, char* argv[]) { - - if ( argc!=3 ) - { - std::cout << std::endl; - std::cout << "Perform Tensor estimation on an dwi file" << std::endl; - std::cout << std::endl; - std::cout << "usage: " << argv[0] << " " << std::endl; - std::cout << std::endl; - return EXIT_FAILURE; - } - - mitk::BaseData::Pointer baseData = 0; - - try - { - mitk::CoreObjectFactory::GetInstance(); - - RegisterDiffusionCoreObjectFactory(); - RegisterCoreExtObjectFactory() ; - - std::string filename = argv[1]; - const std::string s1="", s2=""; - std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false ); - if( infile.empty() ) + if ( argc!=4 ) { - MITK_INFO << "File could not be read" ; + std::cout << std::endl; + std::cout << "Perform Tensor estimation on an dwi file" << std::endl; + std::cout << std::endl; + std::cout << "usage: " << argv[0] << " " << argv[1] << " " << std::endl; + std::cout << std::endl; + return EXIT_FAILURE; } - baseData = infile.at(0); - } - catch ( itk::ExceptionObject &err) - { - MITK_INFO << "Exception during read: " << err; - } - catch ( std::exception err) - { - MITK_INFO << "Exception during read: " << err.what(); - } - catch ( ... ) - { - MITK_INFO << "Exception during read!"; - } - - - - - - typedef short DiffusionPixelType; - typedef float TensorPixelType; - - typedef mitk::DiffusionImage DiffusionImageType; - - - std::string outfilename = argv[2]; - - DiffusionImageType::Pointer vols; - - try - { - vols = dynamic_cast(baseData.GetPointer()); - } - catch ( itk::ExceptionObject &err) - { - MITK_INFO << "Exception during Image write: " << err; - } - catch ( std::exception err) - { - MITK_INFO << "Exception during Image write: " << err.what(); - } - catch ( ... ) - { - MITK_INFO << "Exception during Image write"; - } - - if(vols.IsNull()) - { - return EXIT_FAILURE; - } - - - - typedef float TTensorPixelType; - - typedef itk::DiffusionTensor3DReconstructionImageFilter< - DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; - TensorReconstructionImageFilterType::Pointer filter = - TensorReconstructionImageFilterType::New(); - filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); - filter->SetBValue(vols->GetB_Value()); - filter->Update(); - - - mitk::TensorImage::Pointer image = mitk::TensorImage::New(); - image->InitializeByItk( filter->GetOutput() ); - image->SetVolume( filter->GetOutput()->GetBufferPointer() ); - - - // Save tensor image - std::cout << "writing tensor image"; - - std::string tensorName = outfilename + ".dti"; - mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); - for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) - { - if ( (*it)->CanWriteBaseDataType(image.GetPointer()) ) { - MITK_INFO << "writing"; - (*it)->SetFileName( tensorName.c_str() ); - (*it)->DoWrite( image.GetPointer() ); + try + { + RegisterDiffusionCoreObjectFactory(); + + MITK_INFO << "Loading image ..."; + const std::string s1="", s2=""; + std::vector infile = BaseDataIO::LoadBaseDataFromFile( argv[2], s1, s2, false ); + DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); + + std::string outfilename = argv[3]; + + typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, float > TensorReconstructionImageFilterType; + TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); + filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); + filter->SetBValue(dwi->GetB_Value()); + filter->Update(); + + // Save tensor image + MITK_INFO << "writing image " << outfilename; + itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); + io->SetFileType( itk::ImageIOBase::Binary ); + io->UseCompressionOn(); + + itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< float >, 3 > >::Pointer writer = itk::ImageFileWriter< itk::Image< itk::DiffusionTensor3D< float >, 3 > >::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(outfilename); + writer->SetImageIO(io); + writer->UseCompressionOn(); + writer->Update(); } - } - - - return EXIT_SUCCESS; + catch ( itk::ExceptionObject &err) + { + MITK_INFO << "Exception: " << err; + } + catch ( std::exception err) + { + MITK_INFO << "Exception: " << err.what(); + } + catch ( ... ) + { + MITK_INFO << "Exception!"; + } + return EXIT_SUCCESS; } RegisterDiffusionCoreMiniApp(TensorReconstruction);