diff --git a/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp b/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp index bc66c7b026..3de4c4e693 100644 --- a/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp +++ b/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp @@ -1,370 +1,368 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2007-12-11 14:46:19 +0100 (Di, 11 Dez 2007) $ Version: $Revision: 13129 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkDicomDiffusionImageHeaderReader.h" #include "mitkGEDicomDiffusionImageHeaderReader.h" #include "mitkPhilipsDicomDiffusionImageHeaderReader.h" #include "mitkSiemensDicomDiffusionImageHeaderReader.h" #include "mitkSiemensMosaicDicomDiffusionImageHeaderReader.h" void InsertUnique( std::vector & vec, float value ) { int n = vec.size(); if (n == 0) { vec.push_back( value ); return; } for (int k = 0; k < n ; k++) { if (vec[k] == value) { return; } } // if we get here, it means value is not in vec. vec.push_back( value ); return; } mitk::DicomDiffusionImageHeaderReader::DicomDiffusionImageHeaderReader() { m_SliceOrderIS = true; m_SingleSeries = true; } mitk::DicomDiffusionImageHeaderReader::~DicomDiffusionImageHeaderReader() { } -mitk::DicomDiffusionImageHeaderReader::SupportedVendors +mitk::DicomDiffusionImageHeaderReader::SupportedVendors mitk::DicomDiffusionImageHeaderReader::GetVendorID() { // adapted from namic-sandbox // DicomToNrrdConverter.cxx m_GdcmIO = ImageIOType::New(); m_GdcmIO->LoadPrivateTagsOn(); m_GdcmIO->SetMaxSizeLoadEntry( 65536 ); m_VolumeReader = VolumeReaderType::New(); m_VolumeReader->SetImageIO( m_GdcmIO ); m_VolumeReader->SetFileNames( m_DicomFilenames ); try { m_VolumeReader->Update(); } catch (itk::ExceptionObject &excp) { std::cerr << "Exception thrown while reading slice" << std::endl; std::cerr << excp << std::endl; return SV_UNKNOWN_VENDOR; } - VolumeReaderType::DictionaryArrayRawPointer inputDict + VolumeReaderType::DictionaryArrayRawPointer inputDict = m_VolumeReader->GetMetaDataDictionaryArray(); std::string vendor; itk::ExposeMetaData ( *(*inputDict)[0], "0008|0070", vendor ); // std::cout << vendor << std::endl; std::string ImageType; itk::ExposeMetaData ( *(*inputDict)[0], "0008|0008", ImageType ); //std::cout << ImageType << std::endl; if ( vendor.find("GE") != std::string::npos ) { // for GE data return SV_GE; } else if( vendor.find("SIEMENS") != std::string::npos ) { if ( ImageType.find("MOSAIC") != std::string::npos ) { // for Siemens MOSAIC return SV_SIEMENS_MOSAIC; } else { // for Siemens SPLIT return SV_SIEMENS; } } else if( vendor.find("PHILIPS") != std::string::npos ) { // for philips data return SV_PHILIPS; } else { // for unrecognized vendors return SV_UNKNOWN_VENDOR; } } // do the work void mitk::DicomDiffusionImageHeaderReader::Update() { // check if there are filenames if(m_DicomFilenames.size()) { m_Output = mitk::DiffusionImageHeaderInformation::New(); m_Output->m_DicomFilenames = m_DicomFilenames; // create correct reader switch(GetVendorID()) { case(SV_GE): { GEDicomDiffusionImageHeaderReader::Pointer reader = GEDicomDiffusionImageHeaderReader::New(); reader->SetSeriesDicomFilenames(this->m_DicomFilenames); reader->SetGdcmIO(this->m_GdcmIO); reader->SetVolumeReader(this->m_VolumeReader); reader->SetOutputPointer(this->m_Output); reader->Update(); this->m_Output = reader->GetOutput(); break; } case(SV_SIEMENS): { SiemensDicomDiffusionImageHeaderReader::Pointer reader = SiemensDicomDiffusionImageHeaderReader::New(); reader->SetSeriesDicomFilenames(this->m_DicomFilenames); reader->SetGdcmIO(this->m_GdcmIO); reader->SetVolumeReader(this->m_VolumeReader); reader->SetOutputPointer(this->m_Output); reader->Update(); this->m_Output = reader->GetOutput(); break; } case(SV_SIEMENS_MOSAIC): { SiemensMosaicDicomDiffusionImageHeaderReader::Pointer reader = SiemensMosaicDicomDiffusionImageHeaderReader::New(); reader->SetSeriesDicomFilenames(this->m_DicomFilenames); reader->SetGdcmIO(this->m_GdcmIO); reader->SetVolumeReader(this->m_VolumeReader); reader->SetOutputPointer(this->m_Output); reader->Update(); this->m_Output = reader->GetOutput(); break; } case(SV_PHILIPS): { PhilipsDicomDiffusionImageHeaderReader::Pointer reader = PhilipsDicomDiffusionImageHeaderReader::New(); reader->SetSeriesDicomFilenames(this->m_DicomFilenames); reader->SetGdcmIO(this->m_GdcmIO); reader->SetVolumeReader(this->m_VolumeReader); reader->SetOutputPointer(this->m_Output); reader->Update(); this->m_Output = reader->GetOutput(); break; } case(SV_UNKNOWN_VENDOR): { std::cerr << "diffusion header reader: unknown vendor" << std::endl; break; } } } } // return output -mitk::DiffusionImageHeaderInformation::Pointer +mitk::DiffusionImageHeaderInformation::Pointer mitk::DicomDiffusionImageHeaderReader::GetOutput() { return m_Output; } void mitk::DicomDiffusionImageHeaderReader::ReadPublicTags() { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } - VolumeReaderType::DictionaryArrayRawPointer inputDict + VolumeReaderType::DictionaryArrayRawPointer inputDict = m_VolumeReader->GetMetaDataDictionaryArray(); // load in all public tags - m_nSlice = inputDict->size(); + m_nSlice = inputDict->size(); std::string tag; tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0008|103e", tag ); this->m_Output->seriesDescription = tag.c_str(); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0020|0011", tag ); this->m_Output->seriesNumber = atoi(tag.c_str()); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0010|0010", tag ); this->m_Output->patientName = tag.c_str(); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0028|0010", tag ); this->m_Output->nRows = atoi( tag.c_str() ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0028|0011", tag ); this->m_Output->nCols = atoi( tag.c_str() ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0028|0030", tag ); sscanf( tag.c_str(), "%f\\%f", &this->m_Output->xRes, &this->m_Output->yRes ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0020|0032", tag ); sscanf( tag.c_str(), "%f\\%f\\%f", &this->m_Output->xOrigin, &this->m_Output->yOrigin, &this->m_Output->zOrigin ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0018|0050", tag ); this->m_Output->sliceThickness = atof( tag.c_str() ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0018|0088", tag ); this->m_Output->sliceSpacing = atof( tag.c_str() ); // figure out how many slices are there in a volume, each unique - // SliceLocation represent one slice + // SliceLocation represent one slice for (int k = 0; k < m_nSlice; k++) { tag.clear(); itk::ExposeMetaData ( *(*inputDict)[k], "0020|1041", tag); float sliceLocation = atof( tag.c_str() ); InsertUnique( m_sliceLocations, sliceLocation ); - } + } // check ImageOrientationPatient and figure out slice direction in // L-P-I (right-handed) system. // In Dicom, the coordinate frame is L-P by default. Look at // http://medical.nema.org/dicom/2007/07_03pu.pdf , page 301 tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0020|0037", tag ); float xRow, yRow, zRow, xCol, yCol, zCol, xSlice, ySlice, zSlice /*, orthoSliceSpacing*/; sscanf( tag.c_str(), "%f\\%f\\%f\\%f\\%f\\%f", &xRow, &yRow, &zRow, &xCol, &yCol, &zCol ); //std::cout << "ImageOrientationPatient (0020:0037): "; //std::cout << xRow << ", " << yRow << ", " << zRow << "; "; //std::cout << xCol << ", " << yCol << ", " << zCol << "\n"; // In Dicom, the measurement frame is L-P by default. Look at // http://medical.nema.org/dicom/2007/07_03pu.pdf , page 301, in // order to make this compatible with Slicer's RAS frame, we // multiply the direction cosines by the negatives of the resolution // (resolution is required by nrrd format). Direction cosine is not // affacted since the resulting frame is still a right-handed frame. xRow = -xRow; yRow = -yRow; xCol = -xCol; yCol = -yCol; // Cross product, this gives I-axis direction xSlice = (yRow*zCol-zRow*yCol)*this->m_Output->sliceSpacing; ySlice = (zRow*xCol-xRow*zCol)*this->m_Output->sliceSpacing; zSlice = (xRow*yCol-yRow*xCol)*this->m_Output->sliceSpacing; xRow *= this->m_Output->xRes; yRow *= this->m_Output->xRes; zRow *= this->m_Output->xRes; xCol *= this->m_Output->yRes; yCol *= this->m_Output->yRes; zCol *= this->m_Output->yRes; this->m_Output->xRow = xRow; this->m_Output->yRow = yRow; this->m_Output->zRow = zRow; this->m_Output->xCol = xCol; this->m_Output->yCol = yCol; this->m_Output->zCol = zCol; this->m_Output->xSlice = xSlice; this->m_Output->ySlice = ySlice; this->m_Output->zSlice = zSlice; try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } } -// nRows, nCols, xRes, yRes, xOrigin,yOrigin, -// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, -// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], +// nRows, nCols, xRes, yRes, xOrigin,yOrigin, +// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, +// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], // vendor,SliceMosaic void mitk::DicomDiffusionImageHeaderReader::ReadPublicTags2() { if (!m_SliceOrderIS) { this->m_Output->xSlice = -this->m_Output->xSlice; this->m_Output->ySlice = -this->m_Output->ySlice; this->m_Output->zSlice = -this->m_Output->zSlice; } //std::cout << "Row: " << this->m_Output->xRow << ", " << this->m_Output->yRow << ", " << this->m_Output->zRow << std::endl; //std::cout << "Col: " << this->m_Output->xCol << ", " << this->m_Output->yCol << ", " << this->m_Output->zCol << std::endl; //std::cout << "Sli: " << this->m_Output->xSlice << ", " << this->m_Output->ySlice << ", " << this->m_Output->zSlice << std::endl; //orthoSliceSpacing = fabs(zSlice); //int nVolume; //float maxBvalue = 0; //int nBaseline = 0; } void mitk::DicomDiffusionImageHeaderReader::TransformGradients() { - // transform gradient directions into RAS frame + // transform gradient directions into RAS frame if ( !m_SliceOrderIS ) { this->m_Output->DiffusionVector[2] = -this->m_Output->DiffusionVector[2]; // I -> S } } diff --git a/Modules/DiffusionImaging/DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp b/Modules/DiffusionImaging/DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp index 3da285bbd1..6cccaad051 100644 --- a/Modules/DiffusionImaging/DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp +++ b/Modules/DiffusionImaging/DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp @@ -1,163 +1,161 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2007-12-11 14:46:19 +0100 (Di, 11 Dez 2007) $ Version: $Revision: 13129 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkGEDicomDiffusionImageHeaderReader.h" #include "gdcmGlobal.h" //#include "gdcmVersion.h" #include "gdcmDict.h" #include "gdcmDicts.h" #include "gdcmDictEntry.h" #include "gdcmDictEntry.h" #include "gdcmDict.h" #include "gdcmFile.h" #include "gdcmSerieHelper.h" //const gdcm::DictEntry GEDictBValue( "0043,1039", gdcm::VR::IS, gdcm::VM::VM1, "B Value of diffusion weighting" ); //const gdcm::DictEntry GEDictXGradient( "0019,10bb", gdcm::VR::DS, gdcm::VM::VM1 , "X component of gradient direction" ); //const gdcm::DictEntry GEDictYGradient( "0019,10bc", gdcm::VR::DS, gdcm::VM::VM1 , "Y component of gradient direction" ); //const gdcm::DictEntry GEDictZGradient( "0019,10bd", gdcm::VR::DS, gdcm::VM::VM1 , "Z component of gradient direction" ); mitk::GEDicomDiffusionImageHeaderReader::GEDicomDiffusionImageHeaderReader() { } mitk::GEDicomDiffusionImageHeaderReader::~GEDicomDiffusionImageHeaderReader() { } // do the work void mitk::GEDicomDiffusionImageHeaderReader::Update() { // check if there are filenames if(m_DicomFilenames.size()) { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } // adapted from namic-sandbox // DicomToNrrdConverter.cxx - VolumeReaderType::DictionaryArrayRawPointer inputDict + VolumeReaderType::DictionaryArrayRawPointer inputDict = m_VolumeReader->GetMetaDataDictionaryArray(); ReadPublicTags(); //int mMosaic; // number of raws in each mosaic block; //int nMosaic; // number of columns in each mosaic block float x0, y0, z0; float x1, y1, z1; std::string tag; tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0020|0032", tag ); sscanf( tag.c_str(), "%f\\%f\\%f", &x0, &y0, &z0 ); std::cout << "Slice 0: " << tag << std::endl; tag.clear(); // assume volume interleaving, i.e. the second dicom file stores // the second slice in the same volume as the first dicom file itk::ExposeMetaData ( *(*inputDict)[1], "0020|0032", tag ); sscanf( tag.c_str(), "%f\\%f\\%f", &x1, &y1, &z1 ); std::cout << "Slice 1: " << tag << std::endl; x1 -= x0; y1 -= y0; z1 -= z0; x0 = x1*this->m_Output->xSlice + y1*this->m_Output->ySlice + z1*this->m_Output->zSlice; if (x0 < 0) { m_SliceOrderIS = false; } ReadPublicTags2(); int nSliceInVolume; int nVolume; nSliceInVolume = m_sliceLocations.size(); nVolume = m_nSlice/nSliceInVolume; // assume volume interleaving std::cout << "Number of Slices: " << m_nSlice << std::endl; std::cout << "Number of Volume: " << nVolume << std::endl; std::cout << "Number of Slices in each volume: " << nSliceInVolume << std::endl; for (int k = 0; k < m_nSlice; k += nSliceInVolume) { tag.clear(); bool exist = itk::ExposeMetaData ( *(*inputDict)[k], "0043|1039", tag); float b = atof( tag.c_str() ); this->m_Output->bValue = b; vnl_vector_fixed vect3d; if (!exist || b == 0) { vect3d.fill( 0 ); - this->m_Output->DiffusionVector = vect3d; + this->m_Output->DiffusionVector = vect3d; continue; } vect3d.fill( 0 ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[k], "0019|10bb", tag); vect3d[0] = atof( tag.c_str() ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[k], "0019|10bc", tag); vect3d[1] = atof( tag.c_str() ); tag.clear(); itk::ExposeMetaData ( *(*inputDict)[k], "0019|10bd", tag); vect3d[2] = atof( tag.c_str() ); vect3d.normalize(); - this->m_Output->DiffusionVector = vect3d; + this->m_Output->DiffusionVector = vect3d; } TransformGradients(); try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } } } -//header = new mitk::DWIHeader(nRows, nCols, xRes, yRes, xOrigin,yOrigin, -// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, -// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], +//header = new mitk::DWIHeader(nRows, nCols, xRes, yRes, xOrigin,yOrigin, +// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, +// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], // vendor,SliceMosaic); diff --git a/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp b/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp index 738bf0d461..017e089101 100644 --- a/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp +++ b/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp @@ -1,614 +1,612 @@ /*========================================================================= Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkSiemensDicomDiffusionImageHeaderReader.h" #include "gdcmGlobal.h" //#include #include "gdcmFile.h" #include "gdcmImageReader.h" #include "gdcmDictEntry.h" #include "gdcmDicts.h" #include "gdcmTag.h" mitk::SiemensDicomDiffusionImageHeaderReader::SiemensDicomDiffusionImageHeaderReader() { } mitk::SiemensDicomDiffusionImageHeaderReader::~SiemensDicomDiffusionImageHeaderReader() { } int mitk::SiemensDicomDiffusionImageHeaderReader::ExtractSiemensDiffusionInformation( std::string tagString, std::string nameString, std::vector& valueArray, int startPos ) { std::string::size_type atPosition = tagString.find( nameString, startPos ); if ( atPosition == std::string::npos) { return 0; } else { std::string infoAsString = tagString.substr( atPosition, tagString.size()-atPosition+1 ); const char * infoAsCharPtr = infoAsString.c_str(); int vm = *(infoAsCharPtr+64); std::string vr = infoAsString.substr( 68, 4 ); //int syngodt = *(infoAsCharPtr+72); //int nItems = *(infoAsCharPtr+76); //int localDummy = *(infoAsCharPtr+80); int offset = 84; for (int k = 0; k < vm; k++) { int itemLength = *(infoAsCharPtr+offset+4); int strideSize = static_cast (ceil(static_cast(itemLength)/4) * 4); std::string valueString = infoAsString.substr( offset+16, itemLength ); valueArray.push_back( atof(valueString.c_str()) ); offset += 16+strideSize; } return vm; } } int mitk::SiemensDicomDiffusionImageHeaderReader::ExtractSiemensDiffusionGradientInformation( std::string tagString, std::string nameString, std::vector& valueArray ) { - int nItems = 0; + int nItems = 0; std::string::size_type pos = -1; while(nItems != 3) { - nItems = ExtractSiemensDiffusionInformation( tagString, nameString, valueArray, pos+1 ); + nItems = ExtractSiemensDiffusionInformation( tagString, nameString, valueArray, pos+1 ); pos = tagString.find( nameString, pos+1 ); if ( pos == std::string::npos ) { break; } } return nItems; } // do the work void mitk::SiemensDicomDiffusionImageHeaderReader::Update() { // check if there are filenames if(m_DicomFilenames.size()) { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } // adapted from slicer // DicomToNrrdConverter.cxx - VolumeReaderType::DictionaryArrayRawPointer + VolumeReaderType::DictionaryArrayRawPointer inputDict = m_VolumeReader->GetMetaDataDictionaryArray(); // gdcm::DictEntry SiemensDictBValue( "0019,100c", "B Value of diffusion weighting", gdcm::VR::IS, gdcm::VM::VM1 ); // gdcm::DictEntry SiemensDictDiffusionDirection( "0019,100e", "Diffusion Gradient Direction", gdcm::VR::FD, gdcm::VM::VM3 ); // gdcm::DictEntry SiemensDictDiffusionMatrix( "0019,1027", "Diffusion Matrix", gdcm::VR::FD, gdcm::VM::VM6 ); // gdcm::DictEntry SiemensDictShadowInfo( "0029,1010", "Siemens DWI Info", gdcm::VR::OB, gdcm::VM::VM1 ); // //if(gdcm::Global::GetInstance().GetDicts().GetPublicDict().GetDictEntry(gdcm::Tag(0x0019,0x100c)) == 0) // gdcm::Global::GetInstance().GetDicts().GetPublicDict().AddDictEntry(gdcm::Tag(0x0019,0x100c),SiemensDictBValue); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictDiffusionDirection.GetName()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictDiffusionDirection); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictDiffusionMatrix.GetName()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictDiffusionMatrix); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictShadowInfo.GetName()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictShadowInfo); ReadPublicTags(); //int mMosaic; // number of raws in each mosaic block; //int nMosaic; // number of columns in each mosaic block float x0, y0, z0; float x1, y1, z1; std::string tag; tag.clear(); itk::ExposeMetaData ( *(*inputDict)[0], "0020|0032", tag ); sscanf( tag.c_str(), "%f\\%f\\%f", &x0, &y0, &z0 ); //MITK_INFO << "Slice 0: " << tag << std::endl; tag.clear(); // assume volume interleaving, i.e. the second dicom file stores // the second slice in the same volume as the first dicom file if((*inputDict).size() > 1) { itk::ExposeMetaData ( *(*inputDict)[1], "0020|0032", tag ); sscanf( tag.c_str(), "%f\\%f\\%f", &x1, &y1, &z1 ); //MITK_INFO << "Slice 1: " << tag << std::endl; x1 -= x0; y1 -= y0; z1 -= z0; x0 = x1*this->m_Output->xSlice + y1*this->m_Output->ySlice + z1*this->m_Output->zSlice; if (x0 < 0) { m_SliceOrderIS = false; } } ReadPublicTags2(); int nStride = 1; //MITK_INFO << orthoSliceSpacing << std::endl; this->m_Output->nSliceInVolume = m_sliceLocations.size(); //nVolume = nSlice/this->m_Output->nSliceInVolume; //MITK_INFO << "Number of Slices: " << m_nSlice << std::endl; //MITK_INFO << "Number of Volume: " << nVolume << std::endl; //MITK_INFO << "Number of Slices in each volume: " << this->m_Output->nSliceInVolume << std::endl; nStride = this->m_Output->nSliceInVolume; MITK_INFO << m_DicomFilenames[0] << std::endl; MITK_INFO << "Dims " << this->m_Output->nRows << "x" << this->m_Output->nCols << "x" << this->m_Output->nSliceInVolume << " " << std::endl; for (int k = 0; k < m_nSlice; k += nStride ) { gdcm::ImageReader reader; reader.SetFileName( m_DicomFilenames[k].c_str() ); if( !reader.Read() ) { itkExceptionMacro(<< "Cannot read requested file"); } const gdcm::File &f = reader.GetFile(); const gdcm::DataSet &ds = f.GetDataSet(); // gdcm::DataSet ds = header0->GetDataSet(); gdcm::DataSet::ConstIterator it = ds.Begin(); // Copy of the header->content // copy information stored in 0029,1010 into a string for parsing for(; it != ds.End(); ++it) { const gdcm::DataElement &ref = *it; - if (ref.GetTag() == gdcm::Tag(0x0029,0x1010)) { + if (ref.GetTag() == gdcm::Tag(0x0029,0x1010)) { tag = std::string(ref.GetByteValue()->GetPointer(),ref.GetByteValue()->GetLength()); } } // parse B_value from 0029,1010 tag std::vector valueArray(0); vnl_vector_fixed vect3d; int nItems = ExtractSiemensDiffusionInformation(tag, "B_value", valueArray); if (nItems != 1 || valueArray[0] == 0) // did not find enough information { tag.clear(); MITK_INFO << "Reading diffusion info from 0019|100c and 0019|100e tags" << std::endl; bool success = false; for(it = ds.Begin(); it != ds.End(); ++it) { const gdcm::DataElement &ref = *it; if (ref.GetTag() == gdcm::Tag(0x0019,0x100c)) { tag = std::string(ref.GetByteValue()->GetPointer(),ref.GetByteValue()->GetLength()); this->m_Output->bValue = atof( tag.c_str() ); success = true; } } tag.clear(); if(success) { if(this->m_Output->bValue == 0) { MITK_INFO << "BV: 0 (Baseline image)"; continue; } success = false; for(it = ds.Begin(); it != ds.End(); ++it) { const gdcm::DataElement &ref = *it; if (ref.GetTag() == gdcm::Tag(0x0019,0x100e)) { tag = std::string(ref.GetByteValue()->GetPointer(),ref.GetByteValue()->GetLength()); success = true; } } if(success) { memcpy( &vect3d[0], tag.c_str()+0, 8 ); memcpy( &vect3d[1], tag.c_str()+8, 8 ); memcpy( &vect3d[2], tag.c_str()+16, 8 ); vect3d.normalize(); this->m_Output->DiffusionVector = vect3d; TransformGradients(); MITK_INFO << "BV: " << this->m_Output->bValue; MITK_INFO << " GD: " << this->m_Output->DiffusionVector; continue; } } } - else + else { MITK_INFO << "Reading diffusion info from 0029,1010 tag" << std::endl; this->m_Output->bValue = valueArray[0]; if(this->m_Output->bValue == 0) { MITK_INFO << "BV: 0 (Baseline image)"; continue; } valueArray.resize(0); nItems = ExtractSiemensDiffusionGradientInformation(tag, "DiffusionGradientDirection", valueArray); if (nItems == 3) { vect3d[0] = valueArray[0]; vect3d[1] = valueArray[1]; vect3d[2] = valueArray[2]; vect3d.normalize(); this->m_Output->DiffusionVector = vect3d; TransformGradients(); MITK_INFO << "BV: " << this->m_Output->bValue; MITK_INFO << " GD: " << this->m_Output->DiffusionVector; continue; } } MITK_ERROR << "No diffusion info found, forcing to BASELINE image." << std::endl; this->m_Output->bValue = 0.0; vect3d.fill( 0.0 ); this->m_Output->DiffusionVector = vect3d; } try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } } } -//header = new mitk::DWIHeader(nRows, nCols, xRes, yRes, xOrigin,yOrigin, -// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, -// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], +//header = new mitk::DWIHeader(nRows, nCols, xRes, yRes, xOrigin,yOrigin, +// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, +// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], // vendor,SliceMosaic); // //// do the work //void DicomDiffusionImageHeaderReader::Update() //{ // // // check if there are filenames // if(m_DicomFilenames.IsNotNull() // && m_DicomFilenames->size() > 0) // { // // adapted from namic-sandbox // // DicomToNrrdConverter.cxx // // bool SliceOrderIS = true; // bool SingleSeries = true; // -// VolumeReaderType::DictionaryArrayRawPointer inputDict +// VolumeReaderType::DictionaryArrayRawPointer inputDict // = m_VolumeReader->GetMetaDataDictionaryArray(); // // if ( vendor.find("GE") != std::string::npos ) // { // // for GE data // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(GEDictBValue.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(GEDictBValue); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(GEDictXGradient.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(GEDictXGradient); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(GEDictYGradient.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(GEDictYGradient); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(GEDictZGradient.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(GEDictZGradient); // } // else if( vendor.find("SIEMENS") != std::string::npos ) // { // // for Siemens data // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensMosiacParameters.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensMosiacParameters); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictNMosiac.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictNMosiac); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictBValue.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictBValue); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictDiffusionDirection.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictDiffusionDirection); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictDiffusionMatrix.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictDiffusionMatrix); // if(gdcm::Global::GetDicts()->GetDefaultPubDict()->GetEntry(SiemensDictShadowInfo.GetKey()) == 0) // gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(SiemensDictShadowInfo); // // } // else if( vendor.find("PHILIPS") != std::string::npos ) // { // // for philips data // } // else // { // std::cerr << "Unrecognized vendor.\n" << std::endl; // } // // ReadPublicTags(); // // int mMosaic; // number of raws in each mosaic block; // int nMosaic; // number of columns in each mosaic block // int nSliceInVolume; // // // figure out slice order and mosaic arrangement. // if ( vendor.find("GE") != std::string::npos || // (vendor.find("SIEMENS") != std::string::npos && !SliceMosaic) ) // { // float x0, y0, z0; // float x1, y1, z1; // tag.clear(); // itk::ExposeMetaData ( *(*inputDict)[0], "0020|0032", tag ); // sscanf( tag.c_str(), "%f\\%f\\%f", &x0, &y0, &z0 ); // MITK_INFO << "Slice 0: " << tag << std::endl; // tag.clear(); // // // assume volume interleaving, i.e. the second dicom file stores // // the second slice in the same volume as the first dicom file // itk::ExposeMetaData ( *(*inputDict)[1], "0020|0032", tag ); // sscanf( tag.c_str(), "%f\\%f\\%f", &x1, &y1, &z1 ); // MITK_INFO << "Slice 1: " << tag << std::endl; // x1 -= x0; y1 -= y0; z1 -= z0; // x0 = x1*xSlice + y1*ySlice + z1*zSlice; // if (x0 < 0) // { // SliceOrderIS = false; // } // } // else if ( vendor.find("SIEMENS") != std::string::npos && SliceMosaic ) // { // MITK_INFO << "Siemens SliceMosaic......" << std::endl; // // SliceOrderIS = false; // // // for siemens mosaic image, figure out mosaic slice order from 0029|1010 // tag.clear(); // gdcm::File *header0 = new gdcm::File; // gdcm::BinEntry* binEntry; // // header0->SetMaxSizeLoadEntry(65536); // header0->SetFileName( filenames[0] ); // header0->SetLoadMode( gdcm::LD_ALL ); // header0->Load(); // // // copy information stored in 0029,1010 into a string for parsing // gdcm::DocEntry* docEntry = header0->GetFirstEntry(); // while(docEntry) // { // if ( docEntry->GetKey() == "0029|1010" ) // { // binEntry = dynamic_cast ( docEntry ); // int binLength = binEntry->GetFullLength(); // tag.resize( binLength ); // uint8_t * tagString = binEntry->GetBinArea(); // // for (int k = 0; k < binLength; k++) // { // tag[k] = *(tagString+k); // } // break; // } // docEntry = header0->GetNextEntry(); // } // // // parse SliceNormalVector from 0029,1010 tag // std::vector valueArray(0); // int nItems = ExtractSiemensDiffusionInformation(tag, "SliceNormalVector", valueArray); // if (nItems != 3) // did not find enough information // { // MITK_INFO << "Warning: Cannot find complete information on SliceNormalVector in 0029|1010\n"; // MITK_INFO << " Slice order may be wrong.\n"; // } // else if (valueArray[2] > 0) // { // SliceOrderIS = true; // } // // // parse NumberOfImagesInMosaic from 0029,1010 tag // valueArray.resize(0); // nItems = ExtractSiemensDiffusionInformation(tag, "NumberOfImagesInMosaic", valueArray); // if (nItems == 0) // did not find enough information // { // MITK_INFO << "Warning: Cannot find complete information on NumberOfImagesInMosaic in 0029|1010\n"; // MITK_INFO << " Resulting image may contain empty slices.\n"; // } -// else +// else // { // nSliceInVolume = static_cast(valueArray[0]); // mMosaic = static_cast (ceil(sqrt(valueArray[0]))); // nMosaic = mMosaic; // } // MITK_INFO << "Mosaic in " << mMosaic << " X " << nMosaic << " blocks (total number of blocks = " << valueArray[0] << ").\n"; // } // else // { // } // // ReadPublicTags2(); // // //////////////////////////////////////////////////////////// // // vendor dependent tags. // // read in gradient vectors and determin nBaseline and nMeasurement // // if ( vendor.find("GE") != std::string::npos ) // { // nSliceInVolume = sliceLocations.size(); // nVolume = nSlice/nSliceInVolume; // // // assume volume interleaving // MITK_INFO << "Number of Slices: " << nSlice << std::endl; // MITK_INFO << "Number of Volume: " << nVolume << std::endl; // MITK_INFO << "Number of Slices in each volume: " << nSliceInVolume << std::endl; // // for (int k = 0; k < nSlice; k += nSliceInVolume) // { // tag.clear(); // bool exist = itk::ExposeMetaData ( *(*inputDict)[k], "0043|1039", tag); // float b = atof( tag.c_str() ); // bValues.push_back(b); // // vnl_vector_fixed vect3d; // if (!exist || b == 0) // { // vect3d.fill( 0 ); -// DiffusionVectors.push_back(vect3d); -// DiffusionVectorsOrig.push_back(vect3d); +// DiffusionVectors.push_back(vect3d); +// DiffusionVectorsOrig.push_back(vect3d); // continue; // } // // vect3d.fill( 0 ); // tag.clear(); // itk::ExposeMetaData ( *(*inputDict)[k], "0019|10bb", tag); // vect3d[0] = atof( tag.c_str() ); // // tag.clear(); // itk::ExposeMetaData ( *(*inputDict)[k], "0019|10bc", tag); // vect3d[1] = atof( tag.c_str() ); // // tag.clear(); // itk::ExposeMetaData ( *(*inputDict)[k], "0019|10bd", tag); // vect3d[2] = atof( tag.c_str() ); // -// DiffusionVectorsOrig.push_back(vect3d); +// DiffusionVectorsOrig.push_back(vect3d); // vect3d.normalize(); -// DiffusionVectors.push_back(vect3d); +// DiffusionVectors.push_back(vect3d); // } // } // else if ( vendor.find("SIEMENS") != std::string::npos ) // { // // int nStride = 1; // if ( !SliceMosaic ) // { // MITK_INFO << orthoSliceSpacing << std::endl; // nSliceInVolume = sliceLocations.size(); // nVolume = nSlice/nSliceInVolume; // MITK_INFO << "Number of Slices: " << nSlice << std::endl; // MITK_INFO << "Number of Volume: " << nVolume << std::endl; // MITK_INFO << "Number of Slices in each volume: " << nSliceInVolume << std::endl; // nStride = nSliceInVolume; // } // else // { // MITK_INFO << "Data in Siemens Mosaic Format\n"; // nVolume = nSlice; // MITK_INFO << "Number of Volume: " << nVolume << std::endl; // MITK_INFO << "Number of Slices in each volume: " << nSliceInVolume << std::endl; // nStride = 1; // } // // // for (int k = 0; k < nSlice; k += nStride ) // { // // gdcm::File *header0 = new gdcm::File; // gdcm::BinEntry* binEntry; // // header0->SetMaxSizeLoadEntry(65536); // header0->SetFileName( filenames[k] ); // header0->SetLoadMode( gdcm::LD_ALL ); // header0->Load(); // // // copy information stored in 0029,1010 into a string for parsing // gdcm::DocEntry* docEntry = header0->GetFirstEntry(); // while(docEntry) // { // if ( docEntry->GetKey() == "0029|1010" ) // { // binEntry = dynamic_cast ( docEntry ); // int binLength = binEntry->GetFullLength(); // tag.resize( binLength ); // uint8_t * tagString = binEntry->GetBinArea(); // // for (int n = 0; n < binLength; n++) // { // tag[n] = *(tagString+n); // } // break; // } // docEntry = header0->GetNextEntry(); // } // // // parse B_value from 0029,1010 tag // std::vector valueArray(0); // vnl_vector_fixed vect3d; // int nItems = ExtractSiemensDiffusionInformation(tag, "B_value", valueArray); // if (nItems != 1 || valueArray[0] == 0) // did not find enough information // { // MITK_INFO << "Warning: Cannot find complete information on B_value in 0029|1010\n"; // bValues.push_back( 0.0 ); // vect3d.fill( 0.0 ); -// DiffusionVectors.push_back(vect3d); -// DiffusionVectorsOrig.push_back(vect3d); +// DiffusionVectors.push_back(vect3d); +// DiffusionVectorsOrig.push_back(vect3d); // continue; // } -// else +// else // { // bValues.push_back( valueArray[0] ); // } // // // parse DiffusionGradientDirection from 0029,1010 tag // valueArray.resize(0); // nItems = ExtractSiemensDiffusionInformation(tag, "DiffusionGradientDirection", valueArray); // if (nItems != 3) // did not find enough information // { // MITK_INFO << "Warning: Cannot find complete information on DiffusionGradientDirection in 0029|1010\n"; // vect3d.fill( 0 ); -// DiffusionVectors.push_back(vect3d); -// DiffusionVectorsOrig.push_back(vect3d); +// DiffusionVectors.push_back(vect3d); +// DiffusionVectorsOrig.push_back(vect3d); // } -// else +// else // { // vect3d[0] = valueArray[0]; // vect3d[1] = valueArray[1]; // vect3d[2] = valueArray[2]; -// DiffusionVectorsOrig.push_back(vect3d); +// DiffusionVectorsOrig.push_back(vect3d); // vect3d.normalize(); -// DiffusionVectors.push_back(vect3d); +// DiffusionVectors.push_back(vect3d); // int p = bValues.size(); // MITK_INFO << "Image#: " << k << " BV: " << bValues[p-1] << " GD: " << DiffusionVectors[p-1] << std::endl; // } // } // } // else // { // } // // TransformGradients(); // // /////////////////////////////////////////////// // // construct header info // -// header = new mitk::DWIHeader(nRows, nCols, xRes, yRes, xOrigin,yOrigin, -// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, -// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], +// header = new mitk::DWIHeader(nRows, nCols, xRes, yRes, xOrigin,yOrigin, +// zOrigin, sliceThickness, sliceSpacing,nSliceInVolume, xRow, yRow, +// zRow, xCol,yCol, zCol, xSlice, ySlice, zSlice,bValues[0], DiffusionVectors[0], // vendor,SliceMosaic); // // // set m_Output // // } // //} // diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx index f81797a869..c67837951a 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,364 +1,363 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-02-08 11:19:03 +0100 (Fr, 08 Feb 2008) $ Version: $Revision: 11989 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "itkImageRegionIterator.h" #include "itkImageRegionConstIterator.h" #include "mitkImageCast.h" template mitk::DiffusionImage::DiffusionImage() : m_VectorImage(0), m_Directions(0), m_OriginalDirections(0), m_B_Value(-1.0), m_VectorImageAdaptor(0) { MeasurementFrameType mf; for(int i=0; i<3; i++) for(int j=0; j<3; j++) mf[i][j] = 0; for(int i=0; i<3; i++) mf[i][i] = 1; m_MeasurementFrame = mf; } template mitk::DiffusionImage::~DiffusionImage() { } template void mitk::DiffusionImage ::InitializeFromVectorImage() { if(!m_VectorImage || !m_Directions || m_B_Value==-1.0) { MITK_INFO << "DiffusionImage could not be initialized. Set all members first!" << std::endl; return; } // find bzero index int firstZeroIndex = -1; for(GradientDirectionContainerType::ConstIterator it = m_Directions->Begin(); it != m_Directions->End(); ++it) { firstZeroIndex++; GradientDirectionType g = it.Value(); if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) break; } typedef itk::Image ImgType; typename ImgType::Pointer img = ImgType::New(); img->SetSpacing( m_VectorImage->GetSpacing() ); // Set the image spacing img->SetOrigin( m_VectorImage->GetOrigin() ); // Set the image origin img->SetDirection( m_VectorImage->GetDirection() ); // Set the image direction img->SetLargestPossibleRegion( m_VectorImage->GetLargestPossibleRegion()); img->SetBufferedRegion( m_VectorImage->GetLargestPossibleRegion() ); img->Allocate(); int vecLength = m_VectorImage->GetVectorLength(); InitializeByItk( img.GetPointer(), 1, vecLength ); //for(int i=0; i itw (img, img->GetLargestPossibleRegion() ); itw = itw.Begin(); itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); itr = itr.Begin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(firstZeroIndex)); ++itr; ++itw; } // init SetImportVolume(img->GetBufferPointer());//, 0, 0, CopyMemory); //SetVolume( img->GetBufferPointer(), i ); //} m_DisplayIndex = firstZeroIndex; MITK_INFO << "Diffusion-Image successfully initialized."; } template void mitk::DiffusionImage ::SetDisplayIndexForRendering(int displayIndex) { -MITK_INFO << "ALOHA WE JUST RECEIVED A NEW INDEX FROM THE PROPERTIES IN MITK........................" << displayIndex; int index = displayIndex; int vecLength = m_VectorImage->GetVectorLength(); index = index > vecLength-1 ? vecLength-1 : index; if( m_DisplayIndex != index ) { typedef itk::Image ImgType; typename ImgType::Pointer img = ImgType::New(); CastToItkImage(this, img); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw = itw.Begin(); itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); itr = itr.Begin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(index)); ++itr; ++itw; } } m_DisplayIndex = index; } //template //bool mitk::DiffusionImage::RequestedRegionIsOutsideOfTheBufferedRegion() //{ // return false; //} // //template //void mitk::DiffusionImage::SetRequestedRegion(itk::DataObject * /*data*/) //{ //} // //template //void mitk::DiffusionImage::SetRequestedRegionToLargestPossibleRegion() //{ //} // //template //bool mitk::DiffusionImage::VerifyRequestedRegion() //{ // return true; //} //template //void mitk::DiffusionImage::DuplicateIfSingleSlice() //{ // // new image // typename ImageType::Pointer oldImage = m_Image; // m_Image = ImageType::New(); // m_Image->SetSpacing( oldImage->GetSpacing() ); // Set the image spacing // m_Image->SetOrigin( oldImage->GetOrigin() ); // Set the image origin // m_Image->SetDirection( oldImage->GetDirection() ); // Set the image direction // typename ImageType::RegionType region = oldImage->GetLargestPossibleRegion(); // if(region.GetSize(0) == 1) // region.SetSize(0,3); // if(region.GetSize(1) == 1) // region.SetSize(1,3); // if(region.GetSize(2) == 1) // region.SetSize(2,3); // m_Image->SetLargestPossibleRegion( region ); // m_Image->SetVectorLength( m_Directions->size() ); // m_Image->SetBufferedRegion( region ); // m_Image->Allocate(); // // // average image data that corresponds to identical directions // itk::ImageRegionIterator< ImageType > newIt(m_Image, region); // newIt.GoToBegin(); // itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion()); // oldIt.GoToBegin(); // // while(!newIt.IsAtEnd()) // { // newIt.Set(oldIt.Get()); // ++newIt; // ++oldIt; // if(oldIt.IsAtEnd()) // oldIt.GoToBegin(); // } // //} template bool mitk::DiffusionImage::AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision) { GradientDirectionType diff = g1 - g2; return diff.two_norm() < precision; } template void mitk::DiffusionImage::CorrectDKFZBrokenGradientScheme(double precision) { GradientDirectionContainerType::Pointer directionSet = CalcAveragedDirectionSet(precision, m_Directions); if(directionSet->size() < 7) { MITK_INFO << "Too few directions, assuming and correcting DKFZ-bogus sequence details."; double v [7][3] = {{ 0, 0, 0 }, {-0.707057, 0, 0.707057 }, { 0.707057, 0, 0.707057 }, { 0, 0.707057, 0.707057 }, { 0, 0.707057, -0.707057 }, {-0.707057, 0.707057, 0 }, { 0.707057, 0.707057, 0 } }; int i=0; for(GradientDirectionContainerType::Iterator it = m_OriginalDirections->Begin(); it != m_OriginalDirections->End(); ++it) { it.Value().set(v[i++%7]); } ApplyMeasurementFrame(); } } template mitk::DiffusionImage::GradientDirectionContainerType::Pointer mitk::DiffusionImage::CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions) { // save old and construct new direction container GradientDirectionContainerType::Pointer newDirections = GradientDirectionContainerType::New(); // fill new direction container for(GradientDirectionContainerType::ConstIterator gdcitOld = directions->Begin(); gdcitOld != directions->End(); ++gdcitOld) { // already exists? bool found = false; for(GradientDirectionContainerType::ConstIterator gdcitNew = newDirections->Begin(); gdcitNew != newDirections->End(); ++gdcitNew) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { found = true; break; } } // if not found, add it to new container if(!found) { newDirections->push_back(gdcitOld.Value()); } } return newDirections; } template void mitk::DiffusionImage::AverageRedundantGradients(double precision) { GradientDirectionContainerType::Pointer newDirs = CalcAveragedDirectionSet(precision, m_Directions); GradientDirectionContainerType::Pointer newOriginalDirs = CalcAveragedDirectionSet(precision, m_OriginalDirections); // if sizes equal, we do not need to do anything in this function if(m_Directions->size() == newDirs->size() || m_OriginalDirections->size() == newOriginalDirs->size()) return; GradientDirectionContainerType::Pointer oldDirections = m_Directions; GradientDirectionContainerType::Pointer oldOriginalDirections = m_OriginalDirections; m_Directions = newDirs; m_OriginalDirections = newOriginalDirs; // new image typename ImageType::Pointer oldImage = m_VectorImage; m_VectorImage = ImageType::New(); m_VectorImage->SetSpacing( oldImage->GetSpacing() ); // Set the image spacing m_VectorImage->SetOrigin( oldImage->GetOrigin() ); // Set the image origin m_VectorImage->SetDirection( oldImage->GetDirection() ); // Set the image direction m_VectorImage->SetLargestPossibleRegion( oldImage->GetLargestPossibleRegion() ); m_VectorImage->SetVectorLength( m_Directions->size() ); m_VectorImage->SetBufferedRegion( oldImage->GetLargestPossibleRegion() ); m_VectorImage->Allocate(); // average image data that corresponds to identical directions itk::ImageRegionIterator< ImageType > newIt(m_VectorImage, m_VectorImage->GetLargestPossibleRegion()); newIt.GoToBegin(); itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion()); oldIt.GoToBegin(); // initial new value of voxel typename ImageType::PixelType newVec; newVec.SetSize(m_Directions->size()); newVec.AllocateElements(m_Directions->size()); std::vector > dirIndices; for(GradientDirectionContainerType::ConstIterator gdcitNew = m_Directions->Begin(); gdcitNew != m_Directions->End(); ++gdcitNew) { dirIndices.push_back(std::vector(0)); for(GradientDirectionContainerType::ConstIterator gdcitOld = oldDirections->Begin(); gdcitOld != oldDirections->End(); ++gdcitOld) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { dirIndices[gdcitNew.Index()].push_back(gdcitOld.Index()); } } } int ind1 = -1; while(!newIt.IsAtEnd()) { // progress typename ImageType::IndexType ind = newIt.GetIndex(); ind1 = ind.m_Index[2]; // init new vector with zeros newVec.Fill(0.0); // the old voxel value with duplicates typename ImageType::PixelType oldVec = oldIt.Get(); for(unsigned int i=0; i void mitk::DiffusionImage::ApplyMeasurementFrame() { m_Directions = GradientDirectionContainerType::New(); int c = 0; for(GradientDirectionContainerType::ConstIterator gdcit = m_OriginalDirections->Begin(); gdcit != m_OriginalDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); vec = vec.pre_multiply(m_MeasurementFrame); m_Directions->InsertElement(c, vec); c++; } } diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp index c5a76e4d6c..652b858f57 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp @@ -1,547 +1,532 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ Version: $Revision: 18127 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #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 + // 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!"); + itkWarningMacro("Cache regenerated!"); } if (!m_OutputCache) { - itkWarningMacro("Tree cache is empty!"); + itkWarningMacro("Tree 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()) ->SetOriginalDirections(m_OutputCache->GetOriginalDirections()); 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 == "") + 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 { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } - MITK_INFO << "NrrdDiffusionImageReader READING IMAGE INFORMATION"; + 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(); } else if(ext == ".fsl" || ext == ".fslgz") { // create temporary file with correct ending for nifti-io std::string fname3 = m_FileName; 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); } } - MITK_INFO << "NrrdDiffusionImageReader READING HEADER INFORMATION"; 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) { - MITK_INFO << *itKey << " ---> " << metaString; 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 ); m_OriginalDiffusionVectors->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) { - MITK_INFO << *itKey << " ---> " << metaString; readb0 = true; m_B_Value = atof(metaString.c_str()); } else if (itKey->find("measurement frame") != std::string::npos) { - MITK_INFO << *itKey << " ---> " << metaString; 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(readFrame) { - MITK_INFO << "Applying Measurement Frame: (" << xx << "," << xy << "," << xz - << ") (" << yx << "," << yy << "," << yz << ") (" << zx << "," << zy << "," << zz << ")"; - for(int i=0; i vec(3); vec.copy_in(m_DiffusionVectors->ElementAt(i).data_block()); vec = vec.pre_multiply(m_MeasurementFrame); m_DiffusionVectors->ElementAt(i).copy_in(vec.data_block()); } } - MITK_INFO << "Number of gradient images: " - << numberOfGradientImages - << " and Number of reference images: " - << numberOfImages - numberOfGradientImages; - 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"; } - MITK_INFO << "Found " << bval_entries.size() << " b-values and " << bvec_entries.size() << " bvec-entries"; 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); m_DiffusionVectors->InsertElement(i,vec); m_OriginalDiffusionVectors->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; } // This call updates the output information of the associated VesselTreeData outputForCache->SetVectorImage(img); outputForCache->SetB_Value(m_B_Value); outputForCache->SetDirections(m_DiffusionVectors); outputForCache->SetOriginalDirections(m_OriginalDiffusionVectors); 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 { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; 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()); + 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(); + 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) + 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/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp index 53d61fa9c2..d2c81554a3 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp @@ -1,328 +1,326 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-12-10 18:05:13 +0100 (Mi, 10 Dez 2008) $ Version: $Revision: 15922 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef __mitkNrrdDiffusionImageWriter__cpp #define __mitkNrrdDiffusionImageWriter__cpp #include "mitkNrrdDiffusionImageWriter.h" #include "itkMetaDataDictionary.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkNiftiImageIO.h" #include "itkImageFileWriter.h" #include "itksys/SystemTools.hxx" #include #include template mitk::NrrdDiffusionImageWriter::NrrdDiffusionImageWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } template mitk::NrrdDiffusionImageWriter::~NrrdDiffusionImageWriter() {} template void mitk::NrrdDiffusionImageWriter::GenerateData() { m_Success = false; InputType* input = this->GetInput(); if (input == NULL) { itkWarningMacro(<<"Sorry, input to NrrdDiffusionImageWriter is NULL!"); return; } if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } char keybuffer[512]; char valbuffer[512]; std::string tmp; itk::VectorImage::Pointer img = input->GetVectorImage(); img->GetMetaDataDictionary(); //itk::MetaDataDictionary dic = input->GetImage()->GetMetaDataDictionary(); vnl_matrix_fixed measurementFrame = input->GetMeasurementFrame(); if (measurementFrame(0,0) || measurementFrame(0,1) || measurementFrame(0,2) || measurementFrame(1,0) || measurementFrame(1,1) || measurementFrame(1,2) || measurementFrame(2,0) || measurementFrame(2,1) || measurementFrame(2,2)) { sprintf( valbuffer, " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", measurementFrame(0,0), measurementFrame(0,1), measurementFrame(0,2), measurementFrame(1,0), measurementFrame(1,1), measurementFrame(1,2), measurementFrame(2,0), measurementFrame(2,1), measurementFrame(2,2)); itk::EncapsulateMetaData(input->GetVectorImage()->GetMetaDataDictionary(),std::string("measurement frame"),std::string(valbuffer)); } sprintf( valbuffer, "DWMRI"); itk::EncapsulateMetaData(input->GetVectorImage()->GetMetaDataDictionary(),std::string("modality"),std::string(valbuffer)); if(input->GetOriginalDirections()->Size()) { sprintf( valbuffer, "%1f", input->GetB_Value() ); itk::EncapsulateMetaData(input->GetVectorImage()->GetMetaDataDictionary(),std::string("DWMRI_b-value"),std::string(valbuffer)); } for(unsigned int i=0; iGetOriginalDirections()->Size(); i++) { sprintf( keybuffer, "DWMRI_gradient_%04d", i ); /*if(itk::ExposeMetaData(input->GetMetaDataDictionary(), std::string(keybuffer),tmp)) continue;*/ sprintf( valbuffer, "%1f %1f %1f", input->GetOriginalDirections()->ElementAt(i).get(0), input->GetOriginalDirections()->ElementAt(i).get(1), input->GetOriginalDirections()->ElementAt(i).get(2)); itk::EncapsulateMetaData(input->GetVectorImage()->GetMetaDataDictionary(),std::string(keybuffer),std::string(valbuffer)); } typedef itk::VectorImage ImageType; std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".hdwi" || ext == ".dwi") { itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); //io->SetNrrdVectorType( nrrdKindList ); io->SetFileType( itk::ImageIOBase::Binary ); io->UseCompressionOn(); typedef itk::ImageFileWriter WriterType; typename WriterType::Pointer nrrdWriter = WriterType::New(); nrrdWriter->UseInputMetaDataDictionaryOn(); nrrdWriter->SetInput( input->GetVectorImage() ); nrrdWriter->SetImageIO(io); nrrdWriter->SetFileName(m_FileName); nrrdWriter->UseCompressionOn(); nrrdWriter->SetImageIO(io); try { nrrdWriter->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; throw; } } else if (ext == ".fsl" || ext == ".fslgz") { MITK_INFO << "Writing Nifti-Image for FSL"; typename ImageType::Pointer vecimg = input->GetVectorImage(); typedef itk::Image ImageType4D; typename ImageType4D::Pointer img4 = ImageType4D::New(); typename ImageType::SpacingType spacing = vecimg->GetSpacing(); typename ImageType4D::SpacingType spacing4; for(int i=0; i<3; i++) spacing4[i] = spacing[i]; spacing4[3] = 1; img4->SetSpacing( spacing4 ); // Set the image spacing typename ImageType::PointType origin = vecimg->GetOrigin(); typename ImageType4D::PointType origin4; for(int i=0; i<3; i++) origin4[i] = origin[i]; origin4[3] = 0; img4->SetOrigin( origin4 ); // Set the image origin typename ImageType::DirectionType direction = vecimg->GetDirection(); typename ImageType4D::DirectionType direction4; for(int i=0; i<3; i++) for(int j=0; j<3; j++) direction4[i][j] = direction[i][j]; for(int i=0; i<4; i++) direction4[i][3] = 0; for(int i=0; i<4; i++) direction4[3][i] = 0; direction4[3][3] = 1; img4->SetDirection( direction4 ); // Set the image direction typename ImageType::RegionType region = vecimg->GetLargestPossibleRegion(); typename ImageType4D::RegionType region4; typename ImageType::RegionType::SizeType size = region.GetSize(); typename ImageType4D::RegionType::SizeType size4; for(int i=0; i<3; i++) size4[i] = size[i]; size4[3] = vecimg->GetVectorLength(); typename ImageType::RegionType::IndexType index = region.GetIndex(); typename ImageType4D::RegionType::IndexType index4; for(int i=0; i<3; i++) index4[i] = index[i]; index4[3] = 0; region4.SetSize(size4); region4.SetIndex(index4); img4->SetRegions( region4 ); img4->Allocate(); itk::ImageRegionIterator it (vecimg, vecimg->GetLargestPossibleRegion() ); typedef typename ImageType::PixelType VecPixType; for (it = it.Begin(); !it.IsAtEnd(); ++it) { VecPixType vec = it.Get(); typename ImageType::IndexType currentIndex = it.GetIndex(); for(unsigned int ind=0; indSetPixel(index4, vec[ind]); } } // create copy of file with correct ending for mitk std::string fname3 = m_FileName; std::string newext = ext == ".fsl" ? "nii" : "nii.gz"; std::string::iterator itend = fname3.end(); fname3.replace( itend-3, itend, newext); itk::NiftiImageIO::Pointer io4 = itk::NiftiImageIO::New(); typedef itk::VectorImage ImageType; typedef itk::ImageFileWriter WriterType4; typename WriterType4::Pointer nrrdWriter4 = WriterType4::New(); nrrdWriter4->UseInputMetaDataDictionaryOn(); nrrdWriter4->SetInput( img4 ); nrrdWriter4->SetFileName(fname3); nrrdWriter4->UseCompressionOn(); nrrdWriter4->SetImageIO(io4); try { nrrdWriter4->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; throw; } itksys::SystemTools::CopyAFile(fname3.c_str(), m_FileName.c_str()); if(input->GetDirections()->Size()) { std::ofstream myfile; std::string fname = m_FileName; fname += ".bvals"; myfile.open (fname.c_str()); for(unsigned int i=0; iGetDirections()->Size(); i++) { double twonorm = input->GetDirections()->ElementAt(i).two_norm(); myfile << input->GetB_Value()*twonorm*twonorm << " "; } myfile.close(); std::ofstream myfile2; std::string fname2 = m_FileName; fname2 += ".bvecs"; myfile2.open (fname2.c_str()); for(int j=0; j<3; j++) { for(unsigned int i=0; iGetDirections()->Size(); i++) { myfile2 << input->GetDirections()->ElementAt(i).get(j) << " "; } myfile2 << std::endl; } std::ofstream myfile3; std::string fname4 = m_FileName; fname4 += ".ttk"; myfile3.open (fname4.c_str()); for(unsigned int i=0; iGetDirections()->Size(); i++) { for(int j=0; j<3; j++) { myfile3 << input->GetDirections()->ElementAt(i).get(j) << " "; } myfile3 << std::endl; } } } try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } m_Success = true; } template void mitk::NrrdDiffusionImageWriter::SetInput( InputType* diffVolumes ) { this->ProcessObject::SetNthInput( 0, diffVolumes ); } template mitk::DiffusionImage* mitk::NrrdDiffusionImageWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } template std::vector mitk::NrrdDiffusionImageWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".dwi"); possibleFileExtensions.push_back(".hdwi"); possibleFileExtensions.push_back(".fsl"); possibleFileExtensions.push_back(".fslgz"); return possibleFileExtensions; } #endif //__mitkNrrdDiffusionImageWriter__cpp diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundle.cpp index 0e3690424e..18948c3471 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundle.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundle.cpp @@ -1,1865 +1,1863 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ Version: $Revision: 21975 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkFiberBundle.h" #include "mitkGeometry2D.h" #include "mitkVector.h" #include "mitkPlaneGeometry.h" #include #include "mitkPlanarFigure.h" #include "mitkPlanarCircle.h" #include "mitkPlanarRectangle.h" #include "mitkPlanarPolygon.h" #include "mitkPlanarFigureComposite.h" #include //#include /* statics values to define the position of a fiberTractPoint * related to the plane of a ROI */ const int mitk::FiberBundle::TRACTPOINT_BACKFACE = 0; const int mitk::FiberBundle::TRACTPOINT_ON_PLANE = 1; const int mitk::FiberBundle::TRACTPOINT_FRNTFACE = 2; // implementation of all essential methods from superclass mitk::FiberBundle::FiberBundle() { - MITK_INFO << "init FB"; - m_GroupFiberBundle = FiberGroupType::New(); m_TractContainer = ContainerType::New(); //by default set a standard geometry, usually geometry is set by the user on initializing //a mitkFiberBundle Object mitk::Geometry3D::Pointer fbgeometry = mitk::Geometry3D::New(); fbgeometry->SetIdentity(); this->SetGeometry(fbgeometry); /* for debugging only */ m_debugITKContainer = itkStochTractContainerType::New(); } mitk::FiberBundle::~FiberBundle() { - MITK_INFO << "destructed FB"; + } void mitk::FiberBundle::SetBounds(float* b) { m_boundsFB[0] = b[0]; m_boundsFB[1] = b[1]; m_boundsFB[2] = b[2]; } void mitk::FiberBundle::SetBounds(double* b) { m_boundsFB[0] = b[0]; m_boundsFB[1] = b[1]; m_boundsFB[2] = b[2]; } float* mitk::FiberBundle::GetBounds() { return m_boundsFB; } void mitk::FiberBundle::PushPoint(int fiberIndex, ContainerPointType point) { if( (unsigned)fiberIndex >= m_TractContainer->Size() ) { fiberIndex = m_TractContainer->Size(); ContainerTractType::Pointer tract = ContainerTractType::New(); tract->InsertElement(tract->Size(),point); m_TractContainer->InsertElement(fiberIndex, tract); } else if(fiberIndex>=0) { m_TractContainer->ElementAt(fiberIndex)->InsertElement(m_TractContainer->ElementAt(fiberIndex)->Size(), point); } } void mitk::FiberBundle::PushTract(ContainerTractType::Pointer tract) { m_TractContainer->InsertElement(m_TractContainer->Size(), tract); } mitk::FiberBundle::Pointer mitk::FiberBundle::JoinBundle(mitk::FiberBundle::Pointer bundle) { mitk::FiberBundle::Pointer retval = mitk::FiberBundle::New(); retval->SetGeometry(this->GetGeometry()); retval->SetBounds(this->m_boundsFB); retval->InsertBundle(this); retval->InsertBundle(bundle); return retval; } void mitk::FiberBundle::InsertBundle(mitk::FiberBundle::Pointer bundle) { // INCOMPLETE, SHOULD CHECK FOR DIFFERENT GEOMETRIES IN THIS+BUNDLE // DO INDEX1 -> WORLD -> INDEX2 TRANSFORMATION, IF THEY ARE DIFFERENT. int num = this->GetNumTracts(); FiberGroupType::Pointer groupFiberBundle = bundle->GetGroupFiberBundle(); ChildrenListType *fiberList = groupFiberBundle->GetChildren(); for(ChildrenListType::iterator itLst = fiberList->begin(); itLst != fiberList->end(); ++itLst) { itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << "DTI Tract is NULL!!!!! omg"; continue; } dtiTract->SetId(num++); this->addSingleDTITract(dtiTract); } ContainerType::Pointer container = bundle->GetTractContainer(); for(int i=0; iSize();i++) { this->PushTract(container->ElementAt(i)); } } int mitk::FiberBundle::FindTractByEndpoints(mitk::FiberBundle::DTITubeType::Pointer searchTract) { int searchNrPnts = searchTract->GetNumberOfPoints(); mitk::FiberBundle::DTITubeType::PointListType searchPntList = searchTract->GetPoints(); typedef mitk::FiberBundle::DTITubeType::PointType PointType; DTITubePointType searchPointFirst = searchPntList.at(0); PointType searchPos1 = searchPointFirst.GetPosition(); DTITubePointType searchPointLast = searchPntList.at(searchNrPnts-1); PointType searchPos2 = searchPointLast.GetPosition(); FiberGroupType::Pointer groupFiberBundle = this->GetGroupFiberBundle(); ChildrenListType *fiberList = groupFiberBundle->GetChildren(); for(ChildrenListType::iterator itLst = fiberList->begin(); itLst != fiberList->end(); ++itLst) { itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << "DTI Tract is NULL!!!!! omg"; continue; } mitk::FiberBundle::DTITubeType::PointListType dtiPntList = dtiTract->GetPoints(); int fibrNrPnts = dtiTract->GetNumberOfPoints(); DTITubePointType point_first = dtiPntList.at(0); PointType pos1 = point_first.GetPosition(); DTITubePointType point_last = dtiPntList.at(fibrNrPnts-1); PointType pos2 = point_last.GetPosition(); if( (pos1 == searchPos1 && pos2 == searchPos2) || (pos2 == searchPos1 && pos1 == searchPos2) ) { return dtiTract->GetId(); } } return -1; } mitk::FiberBundle::Pointer mitk::FiberBundle::SubstractBundle(mitk::FiberBundle::Pointer bundle) { mitk::FiberBundle::Pointer retval = mitk::FiberBundle::New(); retval->SetGeometry(this->GetGeometry()); mitk::FiberBundle::Pointer largeBundle = bundle; mitk::FiberBundle::Pointer smallBundle = this; MITK_INFO << "large children " << largeBundle->GetGroupFiberBundle()->GetNumberOfChildren() << "!="<GetNumTracts(); MITK_INFO << "small children " << smallBundle->GetGroupFiberBundle()->GetNumberOfChildren() << "!="<GetNumTracts(); if(this->GetGroupFiberBundle()->GetNumberOfChildren() > bundle->GetGroupFiberBundle()->GetNumberOfChildren()) { MITK_INFO << "this is large (" << this->GetNumTracts() << ">" << bundle->GetNumTracts() << ")"; largeBundle = this; smallBundle = bundle; } ContainerType::Pointer container = largeBundle->GetTractContainer(); int counter = 0; FiberGroupType::Pointer groupFiberBundle = largeBundle->GetGroupFiberBundle(); ChildrenListType *fiberList = groupFiberBundle->GetChildren(); MITK_INFO << "large number children " << groupFiberBundle->GetNumberOfChildren(); for(ChildrenListType::iterator itLst = fiberList->begin(); itLst != fiberList->end(); ++itLst) { itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << "DTI Tract is NULL!!!!! omg"; continue; } int delId = smallBundle->FindTractByEndpoints(dtiTract); if( delId == -1 ) { retval->addSingleDTITract(dtiTract); retval->PushTract(container->ElementAt(counter)); } MITK_INFO << "Counter " << counter++; } return retval; } mitk::FiberBundle::ContainerPointType mitk::FiberBundle::GetPoint(int tractIndex, int pointIndex) { if (tractIndex>=0 && tractIndex=0 && pointIndexGetElement(tractIndex)->GetElement(pointIndex); } return NULL; } int mitk::FiberBundle::GetNumTracts() { return m_TractContainer->Size(); } int mitk::FiberBundle::GetNumPoints(int tractIndex) { if ((unsigned)tractIndex>=0 && (unsigned)tractIndexSize()) { return m_TractContainer->GetElement(tractIndex)->Size(); } else { //added by igi return NULL; } } mitk::FiberBundle::ContainerTractType::Pointer mitk::FiberBundle::GetTract(int tractIndex) { ContainerTractType::Pointer out; if ((unsigned)tractIndex>0 && (unsigned)tractIndexSize()) return m_TractContainer->GetElement(tractIndex); return out; } void mitk::FiberBundle::additkStochTractContainer(itkStochTractContainerType::Pointer itkStochCont) { //for debugging only m_debugITKContainer = itkStochCont; //******************** /* transform itkStochContainer to standardized TractContainer */ for (itkStochTractContainerType::ConstIterator itCont = itkStochCont->Begin(); itCont != itkStochCont->End(); ++itCont) { // get fiber of itkContainer itkStochTractType::Pointer tract = itCont->Value(); itkStochTractType::VertexListType::ConstPointer vertexlist = tract->GetVertexList(); //init a desired containerFiber ContainerTractType::Pointer contTract = ContainerTractType::New(); //get points of fiber for( int j=0; j<(int)vertexlist->Size(); ++j) { // put the point of vertex pointList in itkContainer itkStochTractType::VertexType vertex = vertexlist->GetElement(j); //prepare for dtiPoint ContainerPointType contPnt; contPnt[0] = (float) vertex[0]; contPnt[1] = (float) vertex[1]; contPnt[2] = (float) vertex[2]; contTract->InsertElement(contTract->Size(), contPnt); /* coordinate testing*/ /*if ((float)vertex[0] == contPnt[0]) { MITK_INFO << " [OK] ... X "; }else { MITK_INFO << "[FAIL] ... itkX: " << vertex[0] << " " << contPnt[0] << "\n" ; } if ((float)vertex[1] == contPnt[1]) { MITK_INFO << " [OK] ... Y "; }else { MITK_INFO << "[FAIL] ... itkY: " << vertex[1] << " " << contPnt[1] << "\n" ; } if ((float)vertex[2] == contPnt[2]) { MITK_INFO << " [OK] ... Z " << "\n" ; MITK_INFO << " Values X Y Z: " << contPnt[0] << " " << contPnt[1] << " " << contPnt[2] << "\n"; }else { MITK_INFO << "[FAIL] ... itkZ: " << vertex[2] << " " << contPnt[2] << "\n" ; } */ } // add filled fiber to container m_TractContainer->InsertElement(m_TractContainer->Size(), contTract); } } void mitk::FiberBundle::addTractContainer( ContainerType::Pointer tractContainer ) { m_TractContainer = tractContainer; } void mitk::FiberBundle::initFiberGroup() { /* iterate through the tractContainer and store fibers in DTISpatialObjects */ for (ContainerType::ConstIterator itCont = m_TractContainer->Begin(); itCont != m_TractContainer->End(); ++itCont) { // init new dtiTube and DTITubeType::Pointer dtiTube = DTITubeType::New(); DTITubeType::PointListType list; // get current tract of container ContainerTractType::Pointer contTract = itCont->Value(); // iterate through the number of points ... use iterator, no index-output is needed anyway for(ContainerTractType::Iterator itContTrct = contTract->Begin(); itContTrct != contTract->End(); ++itContTrct) { // init DTITube point DTITubePointType p; ContainerPointType cntp = itContTrct->Value(); p.SetPosition(cntp[0], cntp[1], cntp[2]); p.SetRadius(1); mitk::FiberBundle::fiberPostprocessing_setTensorMatrix( &p ); mitk::FiberBundle::fiberPostprocessing_FA( &p ); //mitk::FiberBundle::fiberPostprocessing_setTensorMatrix( &p ); p.AddField(DTITubePointType::GA, 3); p.SetColor(1,0,0,1); list.push_back(p); } // dtiTubes dtiTube->GetProperty()->SetName("dtiTube"); dtiTube->SetId(itCont->Index()); dtiTube->SetPoints(list); m_GroupFiberBundle->AddSpatialObject(dtiTube); } } void mitk::FiberBundle::fiberPostprocessing_setTensorMatrix(DTITubePointType *dtiP) { float tMatrix[6]; tMatrix[0]=0.0f; tMatrix[1]=1.0f; tMatrix[2]=2.0f; tMatrix[3]=3.0f; tMatrix[4]=4.0f; tMatrix[5]=5.0f; dtiP->SetTensorMatrix(tMatrix); } /* this method is a HowTo method, used for debugging only */ void mitk::FiberBundle::fiberPostprocessing_setPoint(DTITubePointType *dtiP, ContainerPointType vrtx) { /* get values of variable referenced by a ponter double vtxIdx1, vtxIdx2, vtxIdx3; vtxIdx1 = * (double*) &vrtx[0]; vtxIdx2 = * (double*) &vrtx[1]; vtxIdx3 = * (double*) &vrtx[2]; dtiP->SetPosition(vtxIdx1, vtxIdx2, vtxIdx3); */ dtiP->SetPosition(vrtx[0], vrtx[1], vrtx[2]); } void mitk::FiberBundle::fiberPostprocessing_FA(DTITubePointType *dtiP) { debug_PrototypeCounter ++; float valFA; if (debug_PrototypeCounter % 10 < 5) { valFA = 1.0; } else if (debug_PrototypeCounter % 10 == 7) { valFA = 0.3; } else if (debug_PrototypeCounter % 10 == 8) { valFA = 0; } else { valFA = 0.5; } dtiP->AddField(DTITubePointType::FA, valFA); } /* methods for high speed perfoermance dispay or live-monitoring of fiber results */ void mitk::FiberBundle::addContainer4speedDisplay(ContainerType::Pointer speedTractContainer) { } void mitk::FiberBundle::addSingleDTITract(mitk::FiberBundle::DTITubeType::Pointer dtiFbrTrct) { //MITK_INFO << " start addSingleDTITract(): " << m_GroupFiberBundle->GetNumberOfChildren(); m_GroupFiberBundle->AddSpatialObject(dtiFbrTrct); //MITK_INFO << "end addSingleDTITract(): " << m_GroupFiberBundle->GetNumberOfChildren(); } mitk::FiberBundle::FiberGroupType::Pointer mitk::FiberBundle::getGroupFiberBundle() { return m_GroupFiberBundle; } std::vector mitk::FiberBundle::getAllIDsInFiberBundle() { std::vector allFBIds; FiberGroupType::Pointer fiberGroup = this->getGroupFiberBundle(); ChildrenListType *FiberList; FiberList = fiberGroup->GetChildren(); for(ChildrenListType::iterator itLst = FiberList->begin(); itLst != FiberList->end(); ++itLst) { itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << "DTI Tract is NULL!!!!! omg"; continue; } allFBIds.push_back(dtiTract->GetId()); } return allFBIds; } mitk::FiberBundle::Pointer mitk::FiberBundle::extractFibersById(std::vector idSet) { mitk::FiberBundle::Pointer resultingFibers = mitk::FiberBundle::New(); resultingFibers->SetGeometry(this->GetGeometry()); resultingFibers->SetBounds(this->GetBounds()); //get Fiber by ID //get childrenList //compare if current FiberID is in idSet FiberGroupType::Pointer fiberGroup = this->getGroupFiberBundle(); ChildrenListType *FiberList; FiberList = fiberGroup->GetChildren(); MITK_INFO << "writing fibers into datastructure:...."; for (int cg=0; cgPushTract( m_TractContainer->GetElement(trctId) ); } MITK_INFO << "init new fiberBundle..."; resultingFibers->initFiberGroup(); MITK_INFO << "init new fiberBundle [DONE]"; /* for(ChildrenListType::iterator itLst = FiberList->begin(); itLst != FiberList->end(); ++itLst) { itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << "DTI Tract is NULL!!!!! omg"; continue; } // MITK_INFO << "current DTI tract id: " << dtiTract->GetId(); std::vector::iterator retFibers = find(idSet.begin(), idSet.end(), dtiTract->GetId()); if (retFibers != idSet.end()) { //MITK_INFO << "Fiber and ID equal: "; DTITubeType::Pointer copyDTITract = this->copySingleDTITract(dtiTract); MITK_INFO << "fibercontainer before adding new tract" << resultingFibers->GetNumTracts(); resultingFibers->addSingleDTITract(copyDTITract); MITK_INFO << "fibercontainer after adding new tract" << resultingFibers->GetNumTracts(); } //else { // MITK_INFO << "Fiber and ID not ident!"; //} //std::set::iterator retFibers //retFibers = idSet.find(dtiTract->GetId()); } */ return resultingFibers; } /* extract fibers using planar Figures */ //mitk::FiberBundle::Pointer mitk::FiberBundle::extractFibersPF(mitk::PlanarFigure::Pointer pf) std::vector mitk::FiberBundle::extractFibersByPF(mitk::PlanarFigure::Pointer pf, std::vector* smaller_set) { // if incoming pf is a pfc mitk::PlanarFigureComposite::Pointer pfcomp= dynamic_cast(pf.GetPointer()); if (!pfcomp.IsNull()) { //PFC switch (pfcomp->getOperationType()) { case 0: { //AND std::vector childResults = this->extractFibersByPF(pfcomp->getChildAt(0), smaller_set); std::vector tmpChild = childResults; for (int i=1; igetNumberOfChildren(); ++i) { tmpChild = extractFibersByPF(pfcomp->getChildAt(i),&tmpChild); } std::vector AND_Assamblage(tmpChild.begin(), tmpChild.end()); return AND_Assamblage; } case 1: { //OR std::vector childResults = this->extractFibersByPF(pfcomp->getChildAt(0), smaller_set); std::sort(childResults.begin(), childResults.end()); std::vector tmp_container(m_TractContainer->Size()); std::vector::iterator end; for (int i=1; igetNumberOfChildren(); ++i) { std::vector tmpChild = extractFibersByPF(pfcomp->getChildAt(i), smaller_set); sort(tmpChild.begin(), tmpChild.end()); end = std::set_union(childResults.begin(), childResults.end(), tmpChild.begin(), tmpChild.end(), tmp_container.begin() ); childResults.assign(tmp_container.begin(), end); } std::vector OR_Assamblage(childResults.begin(), childResults.end()); return OR_Assamblage; } case 2: { //NOT // FIRST AN AND OPERATION std::vector childResults = this->extractFibersByPF(pfcomp->getChildAt(0), smaller_set); std::sort(childResults.begin(), childResults.end()); std::vector tmpChild = childResults; for (int i=1; igetNumberOfChildren(); ++i) { tmpChild = extractFibersByPF(pfcomp->getChildAt(i),&tmpChild); } std::vector AND_Assamblage(tmpChild.begin(), tmpChild.end()); // get IDs of interesting fibers std::vector interesting_fibers; if(!smaller_set) interesting_fibers = this->getAllIDsInFiberBundle(); else interesting_fibers.assign(smaller_set->begin(), smaller_set->end()); std::sort(interesting_fibers.begin(), interesting_fibers.end()); // all interesting fibers that are not in the AND assemblage std::vector tmp_not(interesting_fibers.size()); std::vector::iterator end; end = std::set_difference(interesting_fibers.begin(), interesting_fibers.end(), AND_Assamblage.begin(), AND_Assamblage.end(), tmp_not.begin() ); std::vector NOT_Assamblage(tmp_not.begin(),end); return NOT_Assamblage; } default: MITK_INFO << "we have an UNDEFINED composition... ERROR" ; break; } } else { mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); mitk::PlanarRectangle::Pointer rectName = mitk::PlanarRectangle::New(); mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); if (pf->GetNameOfClass() == circleName->GetNameOfClass() ) { //MITK_INFO << "We have a circle :-) " ; int cntGaps = 0; //init output FiberBundle mitk::FiberBundle::Pointer FB_Clipped = mitk::FiberBundle::New(); FB_Clipped->SetGeometry(this->GetGeometry()); //get geometry containing the plane which we need to calculate the xrossingPoints mitk::Geometry2D::ConstPointer pfgeometry = pf->GetGeometry2D(); const mitk::PlaneGeometry* planeGeometry = dynamic_cast (pfgeometry.GetPointer()); std::vector XingPoints; std::vector fiberIDs; //################################################################# //############### Find Gap, FrontFace BackFace #################### //################################################################# //get fibercontainer and iterate through the fibers FiberGroupType::Pointer fiberGroup = this->getGroupFiberBundle(); //extract single fibers ChildrenListType *FiberList; FiberList = fiberGroup->GetChildren(); // iterate through each single tract //int cntFibId = -1; for(ChildrenListType::iterator itLst = FiberList->begin(); itLst != FiberList->end(); ++itLst) { //MITK_INFO << "***************** NEW FIBER *********************"; //cntFibId++; itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << " could not cast groupFiberBundle child into dtiTract!... LEVEL 4 ERROR"; continue; } if (smaller_set && std::find(smaller_set->begin(), smaller_set->end(), dtiTract->GetId())==smaller_set->end()) continue; //get list of points int fibrNrPnts = dtiTract->GetNumberOfPoints(); mitk::FiberBundle::DTITubeType::PointListType dtiPntList = dtiTract->GetPoints(); // if one fiber is represented by just one point.... // check if this point is on the plane if (fibrNrPnts <= 0) { MITK_INFO << "HyperERROR in mitkFiberBundle.cpp, No point in fiber....check ur algorithm:"; continue; } /* for finding a gap we need to have knowledge of the previous normal * be aware that there can be more than 1 gaps! */ int prevPntFacing = -999999; int currPntFacing = -999999; //double prevFacing = -99999999099; mitk::Point3D prevFiberPntmm; for(int i=0; iGetGeometry()->IndexToWorld(tmpFiberVec, currentFiberPntmm); double faceing = pfgeometry->SignedDistance(currentFiberPntmm); //MITK_INFO << "FacingOutput: " << faceing; /////////////////////////////////////// if (faceing < 0) { //backface currPntFacing = TRACTPOINT_BACKFACE; //MITK_INFO << "is BACKFACE" << currentFiberPntmm ; } else if (faceing > 0) { //frontface currPntFacing = TRACTPOINT_FRNTFACE; // MITK_INFO << "is FRONTFACE" << currentFiberPntmm ; } else if (faceing == 0) { //onplane currPntFacing = TRACTPOINT_ON_PLANE; //MITK_INFO << "is on PLANE" << currentFiberPntmm ; } //////////////////////////////////////// if (currPntFacing == TRACTPOINT_ON_PLANE) { //strike, //calculate distance //TODO SOURCE THIS CONTROLSTRUCTURE OUT if(isPointInSelection(tmpFiberVec,pf)) { DTITubeType::Pointer copyDTITract = copySingleDTITract(dtiTract); //MITK_INFO << "GroupFB extract after copy of Children: " << fiberGroup->GetNumberOfChildren(); FB_Clipped->addSingleDTITract(copyDTITract); } } else { // Point is BACKFACE or FRONTFACE if (i > 0) { //check if there is a gap between previous and current bool isGap = checkForGap(currPntFacing, prevPntFacing); if (isGap == true) { ++cntGaps; // mitk::Vector3D LineDirection; LineDirection[0] = currentFiberPntmm[0] - prevFiberPntmm[0]; LineDirection[1] = currentFiberPntmm[1] - prevFiberPntmm[1]; LineDirection[2] = currentFiberPntmm[2] - prevFiberPntmm[2]; mitk::Line3D Xingline( prevFiberPntmm, LineDirection ); mitk::Point3D intersectionPoint; Vector3D planeNormal = planeGeometry->GetNormal(); planeNormal.Normalize(); Vector3D lineDirection = Xingline.GetDirection(); lineDirection.Normalize(); double t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { } Vector3D diff; diff = planeGeometry->GetOrigin() - Xingline.GetPoint(); t = ( planeNormal * diff ) / t; intersectionPoint = Xingline.GetPoint() + lineDirection * t; // bool successXing = planeGeometry->IntersectionPoint( Xingline, intersectionPoint ); //mitk::Point3D intersectionPointmm; //planeGeometry->IndexToWorld(intersectionPoint,intersectionPointmm ); // MITK_INFO << "PlaneIntersectionPoint: " << intersectionPoint; if (false) { MITK_INFO << " ERROR CALCULATING INTERSECTION POINT.....should not happen!! "; } //TODO more nice if this if is outside... bool pntInROI = isPointInSelection(intersectionPoint,pf); if(pntInROI) { // MITK_INFO << "POINT ALSO in ROI"; // MITK_INFO << "GroupFB extract before copy new object. of Children: " << fiberGroup->GetNumberOfChildren(); /* dtiTubeSpatial::Pointer can not occur in 2 GroupSpatialObjects, therefore we have to copy the whole dtiTract of current List and add the copy to the desired FiberBundle */ // DTITubeType::Pointer copyDTITract = copySingleDTITract(dtiTract); // MITK_INFO << "GroupFB extract after copy of Children: " << fiberGroup->GetNumberOfChildren(); //FB_Clipped->addSingleDTITract(copyDTITract); //MITK_INFO << "GroupFB extract after adding dtiTract to now FB NR. of Children: " << fiberGroup->GetNumberOfChildren(); //MITK_INFO << "DTI Tract ID: " << dtiTract->GetId(); // MITK_INFO << "size of Vector before pushing: " << fiberIDs.size(); fiberIDs.push_back(dtiTract->GetId()); // MITK_INFO << "size of Vector after pushing: " << fiberIDs.size(); //MITK_INFO << "GroupFB extract after adding dtiTract to now FB NR. of Children: " << fiberGroup->GetNumberOfChildren(); break; } } // MITK_INFO << "---no gap---"; } //endif i>0 } //endif Facing // MITK_INFO << "GroupFB extract end Facing, NR. of Children: " << fiberGroup->GetNumberOfChildren(); //update point flag prevPntFacing = currPntFacing; prevFiberPntmm = currentFiberPntmm; } //end for fiberPoints //MITK_INFO ;<< "GroupFB extract end of single tract, NR. of Children: " << fiberGroup->GetNumberOfChildren(); } //end for fiberTracts // MITK_INFO << "GroupFB extract end of iterating through fiberBundle, NR. of Children: " << fiberGroup->GetNumberOfChildren(); // MITK_INFO << "selected " << fiberIDs.size() << " fibers " << " | counted gaps: " << cntGaps; return fiberIDs; } else if (pf->GetNameOfClass() == rectName->GetNameOfClass() ){ MITK_INFO << "We have a rectangle :-) " ; } else if (pf->GetNameOfClass() == polyName->GetNameOfClass() ) { //MITK_INFO << "We have a polygon :-) " ; int cntGaps = 0; //init output FiberBundle mitk::FiberBundle::Pointer FB_Clipped = mitk::FiberBundle::New(); FB_Clipped->SetGeometry(this->GetGeometry()); //get geometry containing the plane which we need to calculate the xrossingPoints mitk::Geometry2D::ConstPointer pfgeometry = pf->GetGeometry2D(); const mitk::PlaneGeometry* planeGeometry = dynamic_cast (pfgeometry.GetPointer()); std::vector XingPoints; std::vector fiberIDs; //################################################################# //############### Find Gap, FrontFace BackFace #################### //################################################################# //get fibercontainer and iterate through the fibers FiberGroupType::Pointer fiberGroup = this->getGroupFiberBundle(); //extract single fibers ChildrenListType *FiberList; FiberList = fiberGroup->GetChildren(); int nrFibrs = fiberGroup->GetNumberOfChildren(); // iterate through each single tract int cntFib = 1; for(ChildrenListType::iterator itLst = FiberList->begin(); itLst != FiberList->end(); ++itLst) { if (cntFib % 500 == 0) { MITK_INFO << "================\n prosessed fibers: " << cntFib << " of " << nrFibrs;; } ++cntFib; //MITK_INFO << "***************** NEW FIBER *********************"; //cntFibId++; itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << " could not cast groupFiberBundle child into dtiTract!... LEVEL 4 ERROR"; continue; } if (smaller_set && std::find(smaller_set->begin(), smaller_set->end(), dtiTract->GetId())==smaller_set->end()) continue; //get list of points int fibrNrPnts = dtiTract->GetNumberOfPoints(); mitk::FiberBundle::DTITubeType::PointListType dtiPntList = dtiTract->GetPoints(); // if one fiber is represented by just one point.... // check if this point is on the plane if (fibrNrPnts <= 0) { MITK_INFO << "HyperERROR in mitkFiberBundle.cpp, No point in fiber....check ur algorithm:"; continue; } /* for finding a gap we need to have knowledge of the previous normal * be aware that there can be more than 1 gaps! */ int prevPntFacing = -999999; int currPntFacing = -999999; //double prevFacing = -99999999099; mitk::Point3D prevFiberPntmm; for(int i=0; iGetGeometry()->IndexToWorld(tmpFiberVec, currentFiberPntmm); double faceing = pfgeometry->SignedDistance(currentFiberPntmm); //MITK_INFO << "FacingOutput: " << faceing; /////////////////////////////////////// if (faceing < 0) { //backface currPntFacing = TRACTPOINT_BACKFACE; //MITK_INFO << "is BACKFACE" << currentFiberPntmm ; } else if (faceing > 0) { //frontface currPntFacing = TRACTPOINT_FRNTFACE; // MITK_INFO << "is FRONTFACE" << currentFiberPntmm ; } else if (faceing == 0) { //onplane currPntFacing = TRACTPOINT_ON_PLANE; //MITK_INFO << "is on PLANE" << currentFiberPntmm ; } //////////////////////////////////////// if (currPntFacing == TRACTPOINT_ON_PLANE) { //strike, //calculate distance //TODO SOURCE THIS CONTROLSTRUCTURE OUT if(isPointInSelection(tmpFiberVec,pf)) { DTITubeType::Pointer copyDTITract = copySingleDTITract(dtiTract); //MITK_INFO << "GroupFB extract after copy of Children: " << fiberGroup->GetNumberOfChildren(); FB_Clipped->addSingleDTITract(copyDTITract); } } else { // Point is BACKFACE or FRONTFACE if (i > 0) { //check if there is a gap between previous and current bool isGap = checkForGap(currPntFacing, prevPntFacing); if (isGap == true) { ++cntGaps; // mitk::Vector3D LineDirection; LineDirection[0] = currentFiberPntmm[0] - prevFiberPntmm[0]; LineDirection[1] = currentFiberPntmm[1] - prevFiberPntmm[1]; LineDirection[2] = currentFiberPntmm[2] - prevFiberPntmm[2]; mitk::Line3D Xingline( prevFiberPntmm, LineDirection ); mitk::Point3D intersectionPoint; Vector3D planeNormal = planeGeometry->GetNormal(); planeNormal.Normalize(); Vector3D lineDirection = Xingline.GetDirection(); lineDirection.Normalize(); double t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { } Vector3D diff; diff = planeGeometry->GetOrigin() - Xingline.GetPoint(); t = ( planeNormal * diff ) / t; intersectionPoint = Xingline.GetPoint() + lineDirection * t; // bool successXing = planeGeometry->IntersectionPoint( Xingline, intersectionPoint ); //mitk::Point3D intersectionPointmm; //planeGeometry->IndexToWorld(intersectionPoint,intersectionPointmm ); // MITK_INFO << "PlaneIntersectionPoint: " << intersectionPoint; if (false) { MITK_INFO << " ERROR CALCULATING INTERSECTION POINT.....should not happen!! "; } //TODO more nice if this if is outside... bool pntInROI = isPointInSelection(intersectionPoint,pf); if(pntInROI) { // MITK_INFO << "POINT ALSO in ROI"; // MITK_INFO << "GroupFB extract before copy new object. of Children: " << fiberGroup->GetNumberOfChildren(); /* dtiTubeSpatial::Pointer can not occur in 2 GroupSpatialObjects, therefore we have to copy the whole dtiTract of current List and add the copy to the desired FiberBundle */ // DTITubeType::Pointer copyDTITract = copySingleDTITract(dtiTract); // MITK_INFO << "GroupFB extract after copy of Children: " << fiberGroup->GetNumberOfChildren(); //FB_Clipped->addSingleDTITract(copyDTITract); //MITK_INFO << "GroupFB extract after adding dtiTract to now FB NR. of Children: " << fiberGroup->GetNumberOfChildren(); //MITK_INFO << "DTI Tract ID: " << dtiTract->GetId(); // MITK_INFO << "size of Vector before pushing: " << fiberIDs.size(); fiberIDs.push_back(dtiTract->GetId()); // MITK_INFO << "size of Vector after pushing: " << fiberIDs.size(); //MITK_INFO << "GroupFB extract after adding dtiTract to now FB NR. of Children: " << fiberGroup->GetNumberOfChildren(); break; } } // MITK_INFO << "---no gap---"; } //endif i>0 } //endif Facing // MITK_INFO << "GroupFB extract end Facing, NR. of Children: " << fiberGroup->GetNumberOfChildren(); //update point flag prevPntFacing = currPntFacing; prevFiberPntmm = currentFiberPntmm; } //end for fiberPoints //MITK_INFO ;<< "GroupFB extract end of single tract, NR. of Children: " << fiberGroup->GetNumberOfChildren(); } //end for fiberTracts // MITK_INFO << "GroupFB extract end of iterating through fiberBundle, NR. of Children: " << fiberGroup->GetNumberOfChildren(); // MITK_INFO << "selected " << fiberIDs.size() << " fibers " << " | counted gaps: " << cntGaps; return fiberIDs; } else { MITK_INFO << "[ERROR] unhandled PlanarFigureType found!"; } } } mitk::FiberBundle::DTITubeType::Pointer mitk::FiberBundle::copySingleDTITract(mitk::FiberBundle::DTITubeType::Pointer currentDtiFiber) { DTITubeType::Pointer newCopyDtiFiber = DTITubeType::New(); //*** get ID newCopyDtiFiber->SetId( currentDtiFiber->GetId() ); //*** get points DTITubeType::PointListType list; //get list of dtiPoints int fibrNrPnts = currentDtiFiber->GetNumberOfPoints(); mitk::FiberBundle::DTITubeType::PointListType dtiPntList = currentDtiFiber->GetPoints(); for(int i=0; iSetPoints(list); std::string dtiname = currentDtiFiber->GetProperty()->GetName(); newCopyDtiFiber->GetProperty()->SetName(dtiname.c_str()); return newCopyDtiFiber; } /* This Method checks * To recognize a gap between the 2 points, the previous must differ from the current * but if the previous is on the plane and the current one not, then mark this situation as no gap * if both, previous and current are on plane, you'll never jump into this Method */ bool mitk::FiberBundle::checkForGap(int crntPntFacing, int prevPntFacing) { bool isGap = false; if (prevPntFacing != TRACTPOINT_ON_PLANE && crntPntFacing != prevPntFacing) { isGap = true; } else if (prevPntFacing == TRACTPOINT_ON_PLANE && crntPntFacing == TRACTPOINT_ON_PLANE) { MITK_INFO << "###########################################"; MITK_INFO << "$$$ HYPER ERROR, LOGIC MALFUNCTION!! previous point and current point in a fiber are recognized as a potential GAP! THIS IS NOT ALLOWED $$$"; MITK_INFO << "###########################################"; } return isGap; } mitk::Point3D mitk::FiberBundle::calculateCrossingPoint(mitk::Point3D pntFrnt, mitk::Point3D pntbck, mitk::PlanarFigure::Pointer pf) { mitk::Point3D pntXing; //################################################################# //######### Calculate intersection of plane and fiber ############# //################################################################# // transform in space :-) // y = k*x +d // k = (y1 - y0)/(x1 - x0) // d = ok // z adoption, take xy ratio to plane intersection and adopt it to z coordinate // z_intersx = (x1 - intersX)/(x1 - x0) * (z1 - z0) + z1 double y; double k; //slope double d; double x0 = pntFrnt[0]; double y0 = pntFrnt[1]; double z0 = pntFrnt[2]; double x1 = pntbck[0]; double y1 = pntbck[1]; double z1 = pntbck[2]; k = (y1 - y0) / (x1 - x0); // if slope == 0 then leave y as it is, change x // if slope == 1 then leave x as it is, change y // if z of p0 and p1 is the same, just take z /* mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); mitk::PlanarRectangle::Pointer rectName = mitk::PlanarRectangle::New(); mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); if (pf->GetNameOfClass() == circleName->GetNameOfClass() ) { MITK_INFO << "We have a circle :-) " ; //controlpoint 1 is middlepoint //controlpoint 2 is radiuspoint mitk::Vector3D V1w = pf->GetWorldControlPoint(0); //centerPoint mitk::Vector3D V2w = pf->GetWorldControlPoint(1); //radiusPoint mitk::Vector3D V1; mitk::Vector3D V2; this->GetGeometry()->WorldToIndex(V1w, V1); this->GetGeometry()->WorldToIndex(V2w, V2); //calculate distance between those 2 and mitk::Point3D distV; distV = V2 - V1; distV[0] = sqrt( pow(distV[0], 2.0) ); distV[1] = sqrt( pow(distV[1], 2.0) ); distV[2] = sqrt( pow(distV[2], 2.0) ); mitk::Point3D distPnt; distPnt = pnt3D - V1; distPnt[0] = sqrt( pow(distPnt[0], 2.0) ); distPnt[1] = sqrt( pow(distPnt[1], 2.0) ); distPnt[2] = sqrt( pow(distPnt[2], 2.0) ); if ( (distPnt <= distV) ) { pntIsInside = true; } return pntIsInside; //compare the result to the distance of all points an a fiber } else if (pf->GetNameOfClass() == rectName->GetNameOfClass() ){ MITK_INFO << "We have a rectangle :-) " ; } else if (pf->GetNameOfClass() == polyName->GetNameOfClass() ) { MITK_INFO << "We have a polygon :-) " ; } */ return pntXing; } bool mitk::FiberBundle::isPointInSelection(mitk::Point3D pnt3D, mitk::PlanarFigure::Pointer pf) { /* TODO needs to be redesigned.... each time because in planarPolygonsection VTK object will be initialized for each point!...PERFORMANCE LACK!!!!!!! */ //calculate distances mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); mitk::PlanarRectangle::Pointer rectName = mitk::PlanarRectangle::New(); mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); if (pf->GetNameOfClass() == circleName->GetNameOfClass() ) { //MITK_INFO << "We have a circle :-) " ; //controlpoint 1 is middlepoint //controlpoint 2 is radiuspoint mitk::Point3D V1w = pf->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = pf->GetWorldControlPoint(1); //radiusPoint mitk::Vector3D V1; mitk::Vector3D V2; // this->GetGeometry()->WorldToIndex(V1w, V1); // this->GetGeometry()->WorldToIndex(V2w, V2); //calculate distance between those 2 and double distPF; distPF = sqrt((double) (V2w[0] - V1w[0]) * (V2w[0] - V1w[0]) + (V2w[1] - V1w[1]) * (V2w[1] - V1w[1]) + (V2w[2] - V1w[2]) * (V2w[2] - V1w[2])); double XdistPnt = sqrt((double) (pnt3D[0] - V1w[0]) * (pnt3D[0] - V1w[0]) + (pnt3D[1] - V1w[1]) * (pnt3D[1] - V1w[1]) + (pnt3D[2] - V1w[2]) * (pnt3D[2] - V1w[2])) ; if ( (distPF > XdistPnt) ) { return true; } return false; //compare the result to the distance of all points an a fiber } else if (pf->GetNameOfClass() == rectName->GetNameOfClass() ){ MITK_INFO << "We have a rectangle :-) " ; } else if (pf->GetNameOfClass() == polyName->GetNameOfClass() ) { //create vtkPolygon using controlpoints from planarFigure polygon vtkPolygon* polygonVtk = vtkPolygon::New(); //get the control points from pf and insert them to vtkPolygon unsigned int nrCtrlPnts = pf->GetNumberOfControlPoints(); // MITK_INFO << "We have a polygon with " << nrCtrlPnts << " controlpoints: " ; for (int i=0; iGetPoints()->InsertNextPoint((double)pf->GetWorldControlPoint(i)[0], (double)pf->GetWorldControlPoint(i)[1], (double)pf->GetWorldControlPoint(i)[2] ); } //prepare everything for using pointInPolygon function double n[3]; polygonVtk->ComputeNormal(polygonVtk->GetPoints()->GetNumberOfPoints(), static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), n); double bounds[6]; polygonVtk->GetPoints()->GetBounds(bounds); double checkIn[3] = {pnt3D[0], pnt3D[1], pnt3D[2]}; int isInPolygon = polygonVtk->PointInPolygon(checkIn, polygonVtk->GetPoints()->GetNumberOfPoints() , static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), bounds, n); // MITK_INFO << "======IsPointOnPolygon:========\n" << isInPolygon << "\n ======================== "; polygonVtk->Delete(); if (isInPolygon == 1) { // MITK_INFO << "return true"; return true; } else if (isInPolygon == 0) { // MITK_INFO << "return false"; return false; } else { MITK_INFO << "&*&*&*&*&*&*&*&*&*&*&*&*&*&*&*&*&*& \n YOUR DRAWN POLYGON DOES IS DEGENERATED AND NOT ACCEPTED \n DRAW A NEW ONE!! HAI CAPITO \n &*&*&*&*&*&*&*&*&*&*&*&*&*&*&*&*&*&"; } } } void mitk::FiberBundle::debug_members() { /* Debug Workflow description * Test1: check if both, groupSpatialObject and Tract container have the same amount of fiberTracts * Test2: check if points in tracts contain the same coordinates * next... to be continued */ /* ######### TEST 1 START ######### */ //************************************* MITK_INFO << " ############### \n *** TEST Equal Amount of Fibers *** \n "; unsigned int groupEls = m_GroupFiberBundle->GetNumberOfChildren(); unsigned int contEls = m_TractContainer->Size(); unsigned int itkcontEls = m_debugITKContainer->Size(); if (groupEls == contEls && contEls == itkcontEls) { MITK_INFO << "[OK] .... Test1 passed. # of Fibers: " << groupEls << " \n ******************** "; } else { MITK_INFO << "[FAIL]: Container and FiberGroup do not contain same amount of fibers! \n "; // MITK_INFO << " Container # of Fibers: " << contEls << " | FiberBundle # of Fibers: " << groupEls << "\n"; MITK_INFO << " # of Fibers m_Cont: " << contEls << " | GroupFibers: " << groupEls << " | ITKCont: " << itkcontEls << "\n"; } /* ######### TEST 1 END ######### */ //*********************************** /* ######### TEST 2 START ######### */ //************************************* /* iterate through itkcontainer*/ itkStochTractContainerType::ConstIterator itITKCnt; itITKCnt = m_debugITKContainer->Begin(); /* extract DTIFiberTracts of the GroupFiberBundle Object */ // all smartPointers to fibers stored in in a ChildrenList ChildrenListType * FiberList; FiberList = m_GroupFiberBundle->GetChildren(); /* iterate through container, itkcontainer groupFiberBundle in one iteration step */ ChildrenListType::iterator itLst; //STL "list" itself provides no index output of current iteration step. itLst = FiberList->begin(); ContainerType::ConstIterator vecIter; for ( vecIter = m_TractContainer->Begin(); vecIter != m_TractContainer->End(); vecIter++ ) { unsigned int itIdx = vecIter->Index(); MITK_INFO << "FiberIteration: " << itIdx << "\n" ; //get single tract of container ContainerTractType::Pointer contTract = vecIter->Value(); int contNrPnts = contTract->Size(); //get singel tract of itkContainer itkStochTractType::Pointer itkcontTract = itITKCnt->Value(); itkStochTractType::VertexListType::ConstPointer itkcontVrtx = itkcontTract->GetVertexList(); int itkcontNrPnts = itkcontVrtx->Size(); /* lists output is SpatialObject, we know we have DTITubeSpacialObjects dynamic cast only likes pointers, no smartpointers, so each dsmartpointer has membermethod .GetPointer() */ itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; DTITubeType::Pointer dtiTract = dynamic_cast (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { MITK_INFO << " ############### *** ERRROR *** ############### \n ############### *** ERRROR *** ############### \n ############### *** ERRROR *** ############### \n "; return; } //get points of tract int fibrNrPnts = dtiTract->GetNumberOfPoints(); DTITubeType::PointListType dtiPntList = dtiTract->GetPoints(); MITK_INFO << " ############### \n *** TEST Equal Amount of Points in Tract *** \n "; if (contNrPnts == itkcontNrPnts && itkcontNrPnts == fibrNrPnts) { MITK_INFO << "[OK] .... Test2 passed. # of Points in fiber: " << fibrNrPnts << " \n ******************** "; } else { MITK_INFO << "[FAIL]: Tracts do not contain same amount of points! \n "; MITK_INFO << " # of Points m_Cont: " << contNrPnts << " | GroupFibers: " << fibrNrPnts << " | ITKCont: " << itkcontNrPnts << "\n"; } //use for()loop with index instead of iterator cuz of accessing more elements, std vectors provide no index output for(int ip=0; ipGetElement(ip); //get point from itkStochContainer itkStochTractType::VertexType itkPnt = itkcontVrtx->GetElement(ip); //get point from dtiGroup DTITubePointType tmpDtiPnt = dtiPntList.at(ip); DTITubePointType::PointType dtiPoint = tmpDtiPnt.GetPosition(); if (tmpcontPnt[0] == (float)itkPnt[0] && (float)itkPnt[0] == (float)dtiPoint[0]) { MITK_INFO << "TractCont | ITKCont | DTIGroup X: " << tmpcontPnt[0] << "...TEST [OK] " << "\n"; }else{ MITK_INFO << "TractCont | ITKCont | DTIGroup X: " << tmpcontPnt[0] << " " << itkPnt[0] << " " << dtiPoint[0] << "...TEST ##### [FAIL] \n"; } if (tmpcontPnt[1] == (float)itkPnt[1] && (float)itkPnt[1] == (float)dtiPoint[1]) { MITK_INFO << "TractCont | ITKCont | DTIGroup Y: " << tmpcontPnt[1] << "...TEST [OK] " << "\n"; }else{ MITK_INFO << "TractCont | ITKCont | DTIGroup Y: " << tmpcontPnt[1] << " " << itkPnt[1] << " " << dtiPoint[1] << "\n"; } if (tmpcontPnt[2] == (float)itkPnt[2] && (float)itkPnt[2] == (float)dtiPoint[2]) { MITK_INFO << "TractCont | ITKCont | DTIGroup Z: " << tmpcontPnt[2] << "...TEST [OK] " << "\n"; }else{ MITK_INFO << "TractCont | ITKCont | DTIGroup Z: " << tmpcontPnt[2] << " " << itkPnt[2] << " " << dtiPoint[2] << "\n"; } } ++itITKCnt; ++itLst; } } vtkPolyData* mitk::FiberBundle::GeneratePolydata() { MITK_INFO << "writing polydata"; //extractn single fibers //in the groupFiberBundle all smartPointers to single fibers are stored in in a ChildrenList mitk::FiberBundle::ChildrenListType * FiberList; FiberList = this->m_GroupFiberBundle->GetChildren(); /* ######## FIBER PREPARATION END ######### */ /* ######## VTK FIBER REPRESENTATION ######## */ //create a vtkPoints object and store the all the brainFiber points in it vtkPoints *vtkpoints = vtkPoints::New(); //in vtkcells all polylines are stored, actually all id's of them are stored vtkCellArray *vtkcells = vtkCellArray::New(); //in some cases a fiber includes just 1 point, so put it in here vtkCellArray *vtkVrtxs = vtkCellArray::New(); //colors and alpha value for each single point, RGBA = 4 components vtkUnsignedCharArray *colorsT = vtkUnsignedCharArray::New(); colorsT->SetNumberOfComponents(4); colorsT->SetName("ColorValues"); vtkDoubleArray *faColors = vtkDoubleArray::New(); faColors->SetName("FaColors"); //vtkDoubleArray *tubeRadius = vtkDoubleArray::New(); //tubeRadius->SetName("TubeRadius"); // iterate through FiberList for(mitk::FiberBundle::ChildrenListType::iterator itLst = FiberList->begin(); itLst != FiberList->end(); ++itLst) { //all points are stored in one vtkpoints list, soooooooo that the lines find their point id to start and end we need some kind of helper index who monitors the current ids for a polyline unsigned long pntIdxHelper = vtkpoints->GetNumberOfPoints(); // lists output is SpatialObject, we know we have DTITubeSpacialObjects // dynamic cast only likes pointers, no smartpointers, so each dsmartpointer has membermethod .GetPointer() itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; mitk::FiberBundle::DTITubeType::Pointer dtiTract = dynamic_cast< mitk::FiberBundle::DTITubeType * > (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { return NULL; } //get list of points int fibrNrPnts = dtiTract->GetNumberOfPoints(); mitk::FiberBundle::DTITubeType::PointListType dtiPntList = dtiTract->GetPoints(); //create a new polyline for a dtiTract //smartpointer vtkPolyLine *polyLine = vtkPolyLine::New(); polyLine->GetPointIds()->SetNumberOfIds(fibrNrPnts); unsigned char rgba[4] = {255,255,255,255}; //tubeRadius->SetNumberOfTuples(fibrNrPnts); //double tbradius = 1;//default value for radius if (fibrNrPnts <= 0) { //this should never occour! but who knows MITK_INFO << "HyperERROR in fiberBundleMapper3D.cpp ...no point in fiberBundle!!! .. check ur trackingAlgorithm"; continue; } for (int i=0; iGetGeometry()->IndexToWorld(indexPnt, worldPnt); double worldFbrPnt[3] = {worldPnt[0], worldPnt[1], worldPnt[2]}; vtkpoints->InsertNextPoint(worldFbrPnt); // tubeRadius->SetTuple1(i,tbradius); //tuple with 1 argument if (fibrNrPnts == 1) { // if there ist just 1 point in a fiber...wtf, but however represent it as a point vtkVertex *vrtx = vtkVertex::New(); vrtx->GetPointIds()->SetNumberOfIds(1); vrtx->GetPointIds()->SetId(i,i+pntIdxHelper); colorsT->InsertNextTupleValue(rgba); vtkVrtxs->InsertNextCell(vrtx); } else { polyLine->GetPointIds()->SetId(i,i+pntIdxHelper); //colorcoding orientation based if (i0) { //nimm nur diff1 mitk::FiberBundle::DTITubePointType nxttmpFiberPntLst = dtiPntList.at(i+1); mitk::FiberBundle::DTITubePointType::PointType nxttmpFiberPnt = nxttmpFiberPntLst.GetPosition(); //double nxttmpvtkPnt[3] = {nxttmpFiberPnt[0], nxttmpFiberPnt[1], nxttmpFiberPnt[2]}; vnl_vector_fixed< double, 3 > tmpPntvtk((double)tmpvtkPnt[0], (double)tmpvtkPnt[1],(double)tmpvtkPnt[2]); vnl_vector_fixed< double, 3 > nxttmpPntvtk(nxttmpFiberPnt[0], nxttmpFiberPnt[1], nxttmpFiberPnt[2]); vnl_vector_fixed< double, 3 > diff; diff = tmpPntvtk - nxttmpPntvtk; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if(i==0) { //explicit handling of startpoint of line //nimm nur diff1 mitk::FiberBundle::DTITubePointType nxttmpFiberPntLst = dtiPntList.at(i+1); mitk::FiberBundle::DTITubePointType::PointType nxttmpFiberPnt = nxttmpFiberPntLst.GetPosition(); //double nxttmpvtkPnt[3] = {nxttmpFiberPnt[0], nxttmpFiberPnt[1], nxttmpFiberPnt[2]}; vnl_vector_fixed< double, 3 > tmpPntvtk((double)tmpvtkPnt[0], (double)tmpvtkPnt[1],(double)tmpvtkPnt[2]); vnl_vector_fixed< double, 3 > nxttmpPntvtk(nxttmpFiberPnt[0], nxttmpFiberPnt[1], nxttmpFiberPnt[2]); vnl_vector_fixed< double, 3 > diff; diff = tmpPntvtk - nxttmpPntvtk; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if(i==fibrNrPnts) { // nimm nur diff2 mitk::FiberBundle::DTITubePointType nxttmpFiberPntLst = dtiPntList.at(i-1); mitk::FiberBundle::DTITubePointType::PointType nxttmpFiberPnt = nxttmpFiberPntLst.GetPosition(); vnl_vector_fixed< double, 3 > tmpPntvtk((double)tmpvtkPnt[0], (double)tmpvtkPnt[1],(double)tmpvtkPnt[2]); vnl_vector_fixed< double, 3 > nxttmpPntvtk(nxttmpFiberPnt[0], nxttmpFiberPnt[1], nxttmpFiberPnt[2]); vnl_vector_fixed< double, 3 > diff; diff = tmpPntvtk - nxttmpPntvtk; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff[2])); rgba[3] = (unsigned char) (255.0); } colorsT->InsertNextTupleValue(rgba); //get FA value float faVal = tmpFiberPntLst.GetField(mitk::FiberBundle::DTITubePointType::FA); //use insertNextValue cuz FA Values are reperesented as a single number (1 Tuple containing 1 parameter) faColors->InsertNextValue((double) faVal); } } vtkcells->InsertNextCell(polyLine); } //vtkcells->InitTraversal(); // Put points and lines together in one polyData structure m_PolyData = vtkPolyData::New(); m_PolyData->SetPoints(vtkpoints); m_PolyData->SetLines(vtkcells); if (vtkVrtxs->GetSize() > 0) { m_PolyData->SetVerts(vtkVrtxs); } m_PolyData->GetPointData()->AddArray(colorsT); m_PolyData->GetPointData()->AddArray(faColors); return m_PolyData; } /* NECESSARY IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { return true; } void mitk::FiberBundle::SetRequestedRegion( itk::DataObject *data ) { } /* TUTORIAL INSERT POINTS / FIBERS to TRACTCONTAINER */ // points and vectors do not need to be initiated but itkVectorContainer /*ContainerPointType pnt1; pnt1[0] = 1.0; pnt1[1] = 2.0; pnt1[2] = 3.0; ContainerPointType pnt2; pnt2[0] = 4.0; pnt2[1] = 5.0; pnt2[2] = 6.0; ContainerTractType tract1; tract1.push_back(pnt1); tract1.push_back(pnt2); ContainerType::Pointer testContainer = ContainerType::New(); unsigned int freeIdx = testContainer->Size(); MITK_INFO << freeIdx << "\n"; testContainer->InsertElement(freeIdx, tract1); //iterate through all fibers stored in container for(ContainerType::ConstIterator itCont = testContainer->Begin(); itCont != testContainer->End(); itCont++) { //get single tract ContainerTractType tmp_fiber = itCont->Value(); // MITK_INFO << tmp_fiber << "\n"; //iterate through all points within a fibertract for(ContainerTractType::iterator itPnt = tmp_fiber.begin(); itPnt != tmp_fiber.end(); ++itPnt) { // get single point with its coordinates ContainerPointType tmp_pntEx = *itPnt; MITK_INFO << tmp_pntEx[0] << "\n"; MITK_INFO << tmp_pntEx[1] << "\n"; MITK_INFO << tmp_pntEx[2] << "\n"; } } ################### DTI FIBERs TUTORIAL ########################### TUTORIAL HOW TO READ POINTS / FIBERS from DTIGroupSpatialObjectContainer assume our dti fibers are stored in m_GroupFiberBundle // all smartPointers to fibers stored in in a ChildrenList ChildrenListType * FiberList; FiberList = m_GroupFiberBundle->GetChildren(); // iterate through container, itkcontainer groupFiberBundle in one iteration step for(ChildrenListType::iterator itLst = FiberList->begin(); itLst != FiberList->end(); ++FiberList) { // lists output is SpatialObject, we know we have DTITubeSpacialObjects // dynamic cast only likes pointers, no smartpointers, so each dsmartpointer has membermethod .GetPointer() itk::SpatialObject<3>::Pointer tmp_fbr; tmp_fbr = *itLst; DTITubeType::Pointer dtiTract = dynamic_cast (tmp_fbr.GetPointer()); if (dtiTract.IsNull()) { return; } //get list of points int fibrNrPnts = dtiTract->GetNumberOfPoints(); DTITubeType::PointListType dtiPntList = dtiTract->GetPoints(); } */ diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleReader.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleReader.cpp index e337546d19..8c5587435f 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleReader.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleReader.cpp @@ -1,389 +1,412 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ Version: $Revision: 18127 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkFiberBundleReader.h" #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundleReader::XML_GEOMETRY = "geometry"; const char* mitk::FiberBundleReader::XML_MATRIX_XX = "xx"; const char* mitk::FiberBundleReader::XML_MATRIX_XY = "xy"; const char* mitk::FiberBundleReader::XML_MATRIX_XZ = "xz"; const char* mitk::FiberBundleReader::XML_MATRIX_YX = "yx"; const char* mitk::FiberBundleReader::XML_MATRIX_YY = "yy"; const char* mitk::FiberBundleReader::XML_MATRIX_YZ = "yz"; const char* mitk::FiberBundleReader::XML_MATRIX_ZX = "zx"; const char* mitk::FiberBundleReader::XML_MATRIX_ZY = "zy"; const char* mitk::FiberBundleReader::XML_MATRIX_ZZ = "zz"; const char* mitk::FiberBundleReader::XML_ORIGIN_X = "origin_x"; const char* mitk::FiberBundleReader::XML_ORIGIN_Y = "origin_y"; const char* mitk::FiberBundleReader::XML_ORIGIN_Z = "origin_z"; const char* mitk::FiberBundleReader::XML_SPACING_X = "spacing_x"; const char* mitk::FiberBundleReader::XML_SPACING_Y = "spacing_y"; const char* mitk::FiberBundleReader::XML_SPACING_Z = "spacing_z"; const char* mitk::FiberBundleReader::XML_SIZE_X = "size_x"; const char* mitk::FiberBundleReader::XML_SIZE_Y = "size_y"; const char* mitk::FiberBundleReader::XML_SIZE_Z = "size_z"; const char* mitk::FiberBundleReader::XML_FIBER_BUNDLE = "fiber_bundle"; const char* mitk::FiberBundleReader::XML_FIBER = "fiber"; const char* mitk::FiberBundleReader::XML_PARTICLE = "particle"; const char* mitk::FiberBundleReader::XML_ID = "id"; const char* mitk::FiberBundleReader::XML_POS_X = "pos_x"; const char* mitk::FiberBundleReader::XML_POS_Y = "pos_y"; const char* mitk::FiberBundleReader::XML_POS_Z = "pos_z"; const char* mitk::FiberBundleReader::XML_NUM_FIBERS = "num_fibers"; const char* mitk::FiberBundleReader::XML_NUM_PARTICLES = "num_particles"; const char* mitk::FiberBundleReader::XML_FIBER_BUNDLE_FILE = "fiber_bundle_file" ; const char* mitk::FiberBundleReader::XML_FILE_VERSION = "file_version" ; const char* mitk::FiberBundleReader::VERSION_STRING = "0.1" ; namespace mitk { void FiberBundleReader ::GenerateData() { MITK_INFO << "Reading fiber bundle"; if ( ( ! m_OutputCache ) ) { Superclass::SetNumberOfRequiredOutputs(0); this->GenerateOutputInformation(); } if (!m_OutputCache) { itkWarningMacro("Tree cache is empty!"); } Superclass::SetNumberOfRequiredOutputs(1); Superclass::SetNthOutput(0, m_OutputCache.GetPointer()); } void FiberBundleReader::GenerateOutputInformation() { m_OutputCache = OutputType::New(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); ext = itksys::SystemTools::LowerCase(ext); + 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; + } + } + if ( m_FileName == "") { } else if (ext == ".fib") { try { TiXmlDocument doc( m_FileName ); doc.LoadFile(); TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); // save this for later hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement(FiberBundleReader::XML_GEOMETRY).Element(); // read geometry mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); // read origin mitk::Point3D origin; double temp = 0; pElem->Attribute(FiberBundleReader::XML_ORIGIN_X, &temp); origin[0] = temp; pElem->Attribute(FiberBundleReader::XML_ORIGIN_Y, &temp); origin[1] = temp; pElem->Attribute(FiberBundleReader::XML_ORIGIN_Z, &temp); origin[2] = temp; geometry->SetOrigin(origin); // read spacing float spacing[3]; pElem->Attribute(FiberBundleReader::XML_SPACING_X, &temp); spacing[0] = temp; pElem->Attribute(FiberBundleReader::XML_SPACING_Y, &temp); spacing[1] = temp; pElem->Attribute(FiberBundleReader::XML_SPACING_Z, &temp); spacing[2] = temp; geometry->SetSpacing(spacing); // read transform vtkMatrix4x4* m = vtkMatrix4x4::New(); pElem->Attribute(FiberBundleReader::XML_MATRIX_XX, &temp); m->SetElement(0,0,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_XY, &temp); m->SetElement(1,0,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_XZ, &temp); m->SetElement(2,0,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_YX, &temp); m->SetElement(0,1,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_YY, &temp); m->SetElement(1,1,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_YZ, &temp); m->SetElement(2,1,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_ZX, &temp); m->SetElement(0,2,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_ZY, &temp); m->SetElement(1,2,temp); pElem->Attribute(FiberBundleReader::XML_MATRIX_ZZ, &temp); m->SetElement(2,2,temp); m->SetElement(0,3,origin[0]); m->SetElement(1,3,origin[1]); m->SetElement(2,3,origin[2]); m->SetElement(3,3,1); geometry->SetIndexToWorldTransformByVtkMatrix(m); // read bounds float bounds[] = {0, 0, 0, 0, 0, 0}; pElem->Attribute(FiberBundleReader::XML_SIZE_X, &temp); bounds[1] = temp; pElem->Attribute(FiberBundleReader::XML_SIZE_Y, &temp); bounds[3] = temp; pElem->Attribute(FiberBundleReader::XML_SIZE_Z, &temp); bounds[5] = temp; geometry->SetFloatBounds(bounds); // read bounds float bounds2[] = {0, 0, 0}; bounds2[0] = bounds[1]; bounds2[1] = bounds[3]; bounds2[2] = bounds[5]; m_OutputCache->SetBounds(bounds2); geometry->SetImageGeometry(true); m_OutputCache->SetGeometry(geometry); // generate tract container ContainerType::Pointer tractContainer = ContainerType::New(); int fiberID = 0; pElem = hRoot.FirstChildElement(FiberBundleReader::XML_FIBER_BUNDLE).FirstChild().Element(); for( pElem; pElem; pElem=pElem->NextSiblingElement()) { TiXmlElement* pElem2 = pElem->FirstChildElement(); ContainerTractType::Pointer tract = ContainerTractType::New(); for( pElem2; pElem2; pElem2=pElem2->NextSiblingElement()) { ContainerPointType point; pElem2->Attribute(FiberBundleReader::XML_POS_X, &temp); point[0] = temp; pElem2->Attribute(FiberBundleReader::XML_POS_Y, &temp); point[1] = temp; pElem2->Attribute(FiberBundleReader::XML_POS_Z, &temp); point[2] = temp; tract->InsertElement(tract->Size(), point); } pElem->Attribute(FiberBundleReader::XML_ID, &fiberID); tractContainer->CreateIndex(fiberID); tractContainer->SetElement(fiberID, tract); } m_OutputCache->addTractContainer(tractContainer); m_OutputCache->initFiberGroup(); MITK_INFO << "Fiber bundle read"; } catch(...) { MITK_INFO << "Could not read file "; } } else if (ext == ".vfib") { // generate tract container ContainerType::Pointer tractContainer = ContainerType::New(); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); ///We create a Generic Reader to test de .vtk/ vtkDataReader *chooser=vtkDataReader::New(); chooser->SetFileName(m_FileName.c_str() ); if( chooser->IsFilePolyData()) { vtkPolyDataReader *reader = vtkPolyDataReader::New(); reader->SetFileName( m_FileName.c_str() ); reader->Update(); if ( reader->GetOutput() != NULL ) { vtkPolyData* output = reader->GetOutput(); output->ComputeBounds(); double bounds[3]; output->GetBounds(bounds); double center[3]; output->GetCenter(center); Point3D origin; origin.SetElement(0, -center[0]); origin.SetElement(1, -center[1]); origin.SetElement(2, -center[2]); - MITK_INFO << origin; mitk::Surface::Pointer surf = mitk::Surface::New(); surf->SetVtkPolyData(output); mitk::Geometry3D* geom = surf->GetGeometry(); //geom->SetOrigin(origin); geom->SetImageGeometry(true); m_OutputCache->SetBounds(bounds); m_OutputCache->SetGeometry(geom); vtkCellArray* cells = output->GetLines(); cells->InitTraversal(); for (int i=0; iGetNumberOfCells(); i++) { ContainerTractType::Pointer tract = ContainerTractType::New(); vtkCell* cell = output->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p); ContainerPointType point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; tract->InsertElement(tract->Size(), point); } tractContainer->InsertElement(i, tract); } } reader->Delete(); } chooser->Delete(); m_OutputCache->addTractContainer(tractContainer); m_OutputCache->initFiberGroup(); MITK_INFO << "Fiber bundle read"; } + + try + { + setlocale(LC_ALL, currLocale.c_str()); + } + catch(...) + { + MITK_INFO << "Could not reset locale " << currLocale; + } } void FiberBundleReader::Update() { this->GenerateData(); } const char* FiberBundleReader ::GetFileName() const { return m_FileName.c_str(); } void FiberBundleReader ::SetFileName(const char* aFileName) { m_FileName = aFileName; } const char* FiberBundleReader ::GetFilePrefix() const { return m_FilePrefix.c_str(); } void FiberBundleReader ::SetFilePrefix(const char* aFilePrefix) { m_FilePrefix = aFilePrefix; } const char* FiberBundleReader ::GetFilePattern() const { return m_FilePattern.c_str(); } void FiberBundleReader ::SetFilePattern(const char* aFilePattern) { m_FilePattern = aFilePattern; } bool FiberBundleReader ::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 == ".fib" || ext == ".vfib") { return true; } return false; } } //namespace MITK diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp index 995d534aff..41f24ed4a5 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp @@ -1,231 +1,260 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-12-10 18:05:13 +0100 (Mi, 10 Dez 2008) $ Version: $Revision: 15922 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkFiberBundleWriter.h" #include const char* mitk::FiberBundleWriter::XML_GEOMETRY = "geometry"; const char* mitk::FiberBundleWriter::XML_MATRIX_XX = "xx"; const char* mitk::FiberBundleWriter::XML_MATRIX_XY = "xy"; const char* mitk::FiberBundleWriter::XML_MATRIX_XZ = "xz"; const char* mitk::FiberBundleWriter::XML_MATRIX_YX = "yx"; const char* mitk::FiberBundleWriter::XML_MATRIX_YY = "yy"; const char* mitk::FiberBundleWriter::XML_MATRIX_YZ = "yz"; const char* mitk::FiberBundleWriter::XML_MATRIX_ZX = "zx"; const char* mitk::FiberBundleWriter::XML_MATRIX_ZY = "zy"; const char* mitk::FiberBundleWriter::XML_MATRIX_ZZ = "zz"; const char* mitk::FiberBundleWriter::XML_ORIGIN_X = "origin_x"; const char* mitk::FiberBundleWriter::XML_ORIGIN_Y = "origin_y"; const char* mitk::FiberBundleWriter::XML_ORIGIN_Z = "origin_z"; const char* mitk::FiberBundleWriter::XML_SPACING_X = "spacing_x"; const char* mitk::FiberBundleWriter::XML_SPACING_Y = "spacing_y"; const char* mitk::FiberBundleWriter::XML_SPACING_Z = "spacing_z"; const char* mitk::FiberBundleWriter::XML_SIZE_X = "size_x"; const char* mitk::FiberBundleWriter::XML_SIZE_Y = "size_y"; const char* mitk::FiberBundleWriter::XML_SIZE_Z = "size_z"; const char* mitk::FiberBundleWriter::XML_FIBER_BUNDLE = "fiber_bundle"; const char* mitk::FiberBundleWriter::XML_FIBER = "fiber"; const char* mitk::FiberBundleWriter::XML_PARTICLE = "particle"; const char* mitk::FiberBundleWriter::XML_ID = "id"; const char* mitk::FiberBundleWriter::XML_POS_X = "pos_x"; const char* mitk::FiberBundleWriter::XML_POS_Y = "pos_y"; const char* mitk::FiberBundleWriter::XML_POS_Z = "pos_z"; const char* mitk::FiberBundleWriter::XML_NUM_FIBERS = "num_fibers"; const char* mitk::FiberBundleWriter::XML_NUM_PARTICLES = "num_particles"; const char* mitk::FiberBundleWriter::XML_FIBER_BUNDLE_FILE = "fiber_bundle_file" ; const char* mitk::FiberBundleWriter::XML_FILE_VERSION = "file_version" ; const char* mitk::FiberBundleWriter::VERSION_STRING = "0.1" ; const char* mitk::FiberBundleWriter::ASCII_FILE = "ascii_file" ; const char* mitk::FiberBundleWriter::FILE_NAME = "file_name" ; mitk::FiberBundleWriter::FiberBundleWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } mitk::FiberBundleWriter::~FiberBundleWriter() {} void mitk::FiberBundleWriter::GenerateData() { - MITK_INFO << "Writing fiber bundle"; - m_Success = false; - InputType* input = this->GetInput(); - if (input == NULL) + try { - itkWarningMacro(<<"Sorry, input to FiberBundleWriter is NULL!"); - return; - } - if ( m_FileName == "" ) - { - itkWarningMacro( << "Sorry, filename has not been set!" ); - return ; - } + MITK_INFO << "Writing fiber bundle"; + m_Success = false; + InputType* input = this->GetInput(); + if (input == NULL) + { + itkWarningMacro(<<"Sorry, input to FiberBundleWriter is NULL!"); + return; + } + if ( m_FileName == "" ) + { + itkWarningMacro( << "Sorry, filename has not been set!" ); + return ; + } - std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); - ext = itksys::SystemTools::LowerCase(ext); + const std::string& locale = "C"; + const std::string& currLocale = setlocale( LC_ALL, NULL ); - if (ext == ".fib") - { - /* direct linked includes of mitkFiberBundle DataStructure */ - typedef mitk::FiberBundle::ContainerPointType ContainerPointType; - typedef mitk::FiberBundle::ContainerTractType ContainerTractType; - typedef mitk::FiberBundle::ContainerType ContainerType; - - ContainerType::Pointer tractContainer = input->GetTractContainer(); - mitk::Geometry3D* geometry = input->GetGeometry(); - - TiXmlDocument documentXML; - TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); // TODO what to write here? encoding? etc.... - documentXML.LinkEndChild( declXML ); - - TiXmlElement* mainXML = new TiXmlElement(mitk::FiberBundleWriter::XML_FIBER_BUNDLE_FILE); - mainXML->SetAttribute(mitk::FiberBundleWriter::XML_FILE_VERSION, mitk::FiberBundleWriter::VERSION_STRING); - documentXML.LinkEndChild(mainXML); - - TiXmlElement* geometryXML = new TiXmlElement(mitk::FiberBundleWriter::XML_GEOMETRY); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_XX, geometry->GetMatrixColumn(0)[0]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_XY, geometry->GetMatrixColumn(0)[1]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_XZ, geometry->GetMatrixColumn(0)[2]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_YX, geometry->GetMatrixColumn(1)[0]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_YY, geometry->GetMatrixColumn(1)[1]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_YZ, geometry->GetMatrixColumn(1)[2]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_ZX, geometry->GetMatrixColumn(2)[0]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_ZY, geometry->GetMatrixColumn(2)[1]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_ZZ, geometry->GetMatrixColumn(2)[2]); - - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_ORIGIN_X, geometry->GetOrigin()[0]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_ORIGIN_Y, geometry->GetOrigin()[1]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_ORIGIN_Z, geometry->GetOrigin()[2]); - - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SPACING_X, geometry->GetSpacing()[0]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SPACING_Y, geometry->GetSpacing()[1]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SPACING_Z, geometry->GetSpacing()[2]); - - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SIZE_X, input->GetBounds()[0]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SIZE_Y, input->GetBounds()[1]); - geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SIZE_Z, input->GetBounds()[2]); - - mainXML->LinkEndChild(geometryXML); - - TiXmlElement* fiberBundleXML = new TiXmlElement(mitk::FiberBundleWriter::XML_FIBER_BUNDLE); - fiberBundleXML->SetAttribute(mitk::FiberBundleWriter::XML_NUM_FIBERS, tractContainer->size()); - int numParticles = 0; - for (int i=0; iSize(); i++) + if ( locale.compare(currLocale)!=0 ) { - ContainerTractType::Pointer tract = tractContainer->GetElement(i); - TiXmlElement* fiberXML = new TiXmlElement(mitk::FiberBundleWriter::XML_FIBER); - fiberXML->SetAttribute(mitk::FiberBundleWriter::XML_ID, i); - fiberXML->SetAttribute(mitk::FiberBundleWriter::XML_NUM_PARTICLES, tract->Size()); - numParticles += tract->Size(); - for (int j=0; jSize(); j++) + try + { + setlocale(LC_ALL, locale.c_str()); + } + catch(...) { - TiXmlElement* particleXML = new TiXmlElement(mitk::FiberBundleWriter::XML_PARTICLE); - ContainerPointType p = tract->GetElement(j); - particleXML->SetAttribute(mitk::FiberBundleWriter::XML_ID, j); - particleXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_POS_X, p[0]); - particleXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_POS_Y, p[1]); - particleXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_POS_Z, p[2]); - fiberXML->LinkEndChild(particleXML); + MITK_INFO << "Could not set locale " << locale; } - fiberBundleXML->LinkEndChild(fiberXML); } - fiberBundleXML->SetAttribute(mitk::FiberBundleWriter::XML_NUM_PARTICLES, numParticles); - mainXML->LinkEndChild(fiberBundleXML); - documentXML.SaveFile( m_FileName ); + std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); + ext = itksys::SystemTools::LowerCase(ext); - m_Success = true; + if (ext == ".fib") + { + /* direct linked includes of mitkFiberBundle DataStructure */ + typedef mitk::FiberBundle::ContainerPointType ContainerPointType; + typedef mitk::FiberBundle::ContainerTractType ContainerTractType; + typedef mitk::FiberBundle::ContainerType ContainerType; + + ContainerType::Pointer tractContainer = input->GetTractContainer(); + mitk::Geometry3D* geometry = input->GetGeometry(); + + TiXmlDocument documentXML; + TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); // TODO what to write here? encoding? etc.... + documentXML.LinkEndChild( declXML ); + + TiXmlElement* mainXML = new TiXmlElement(mitk::FiberBundleWriter::XML_FIBER_BUNDLE_FILE); + mainXML->SetAttribute(mitk::FiberBundleWriter::XML_FILE_VERSION, mitk::FiberBundleWriter::VERSION_STRING); + documentXML.LinkEndChild(mainXML); + + TiXmlElement* geometryXML = new TiXmlElement(mitk::FiberBundleWriter::XML_GEOMETRY); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_XX, geometry->GetMatrixColumn(0)[0]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_XY, geometry->GetMatrixColumn(0)[1]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_XZ, geometry->GetMatrixColumn(0)[2]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_YX, geometry->GetMatrixColumn(1)[0]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_YY, geometry->GetMatrixColumn(1)[1]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_YZ, geometry->GetMatrixColumn(1)[2]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_ZX, geometry->GetMatrixColumn(2)[0]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_ZY, geometry->GetMatrixColumn(2)[1]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_MATRIX_ZZ, geometry->GetMatrixColumn(2)[2]); + + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_ORIGIN_X, geometry->GetOrigin()[0]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_ORIGIN_Y, geometry->GetOrigin()[1]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_ORIGIN_Z, geometry->GetOrigin()[2]); + + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SPACING_X, geometry->GetSpacing()[0]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SPACING_Y, geometry->GetSpacing()[1]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SPACING_Z, geometry->GetSpacing()[2]); + + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SIZE_X, input->GetBounds()[0]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SIZE_Y, input->GetBounds()[1]); + geometryXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_SIZE_Z, input->GetBounds()[2]); + + mainXML->LinkEndChild(geometryXML); + + TiXmlElement* fiberBundleXML = new TiXmlElement(mitk::FiberBundleWriter::XML_FIBER_BUNDLE); + fiberBundleXML->SetAttribute(mitk::FiberBundleWriter::XML_NUM_FIBERS, tractContainer->size()); + int numParticles = 0; + for (int i=0; iSize(); i++) + { + ContainerTractType::Pointer tract = tractContainer->GetElement(i); + TiXmlElement* fiberXML = new TiXmlElement(mitk::FiberBundleWriter::XML_FIBER); + fiberXML->SetAttribute(mitk::FiberBundleWriter::XML_ID, i); + fiberXML->SetAttribute(mitk::FiberBundleWriter::XML_NUM_PARTICLES, tract->Size()); + numParticles += tract->Size(); + for (int j=0; jSize(); j++) + { + TiXmlElement* particleXML = new TiXmlElement(mitk::FiberBundleWriter::XML_PARTICLE); + ContainerPointType p = tract->GetElement(j); + particleXML->SetAttribute(mitk::FiberBundleWriter::XML_ID, j); + particleXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_POS_X, p[0]); + particleXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_POS_Y, p[1]); + particleXML->SetDoubleAttribute(mitk::FiberBundleWriter::XML_POS_Z, p[2]); + fiberXML->LinkEndChild(particleXML); + } + fiberBundleXML->LinkEndChild(fiberXML); + } + fiberBundleXML->SetAttribute(mitk::FiberBundleWriter::XML_NUM_PARTICLES, numParticles); + mainXML->LinkEndChild(fiberBundleXML); - MITK_INFO << "Fiber bundle written"; + documentXML.SaveFile( m_FileName ); - }else if (ext == ".vfib" || ext == ".vtk") { - vtkPolyDataWriter* writer = vtkPolyDataWriter::New(); - writer->SetInput(input->GeneratePolydata()); - writer->SetFileName(m_FileName.c_str()); - writer->SetFileTypeToASCII(); - writer->Write(); + m_Success = true; - m_Success = true; - MITK_INFO << "Fiber bundle written as polydata"; - } + MITK_INFO << "Fiber bundle written"; + + }else if (ext == ".vfib" || ext == ".vtk") { + vtkPolyDataWriter* writer = vtkPolyDataWriter::New(); + writer->SetInput(input->GeneratePolydata()); + writer->SetFileName(m_FileName.c_str()); + writer->SetFileTypeToASCII(); + writer->Write(); + m_Success = true; + MITK_INFO << "Fiber bundle written as polydata"; + } + try + { + setlocale(LC_ALL, currLocale.c_str()); + } + catch(...) + { + MITK_INFO << "Could not reset locale " << currLocale; + } + } + catch(...) + { + throw; + } } void mitk::FiberBundleWriter::SetInputFiberBundle( InputType* diffVolumes ) { this->ProcessObject::SetNthInput( 0, diffVolumes ); } mitk::FiberBundle* mitk::FiberBundleWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } std::vector mitk::FiberBundleWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".fib"); possibleFileExtensions.push_back(".vfib"); possibleFileExtensions.push_back(".vtk"); return possibleFileExtensions; } diff --git a/Modules/DiffusionImaging/IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp b/Modules/DiffusionImaging/IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp index 71ac766641..5292efe764 100644 --- a/Modules/DiffusionImaging/IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp +++ b/Modules/DiffusionImaging/IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp @@ -1,131 +1,119 @@ /* * mitkPlanarFigureComposite.cpp * mitk-all * * Created by HAL9000 on 2/4/11. * Copyright 2011 __MyCompanyName__. All rights reserved. * */ #include "mitkPlanarFigureComposite.h" mitk::PlanarFigureComposite::PlanarFigureComposite() { m_PFVector = CompositionContainer::New(); m_DNVector = DataNodeContainer::New(); } mitk::PlanarFigureComposite::~PlanarFigureComposite() { - + } void mitk::PlanarFigureComposite::addDataNode(mitk::DataNode::Pointer dnode) { - + m_DNVector->InsertElement(m_DNVector->Size(), dnode); - + } void mitk::PlanarFigureComposite::addPlanarFigure(PlanarFigure::Pointer pf) { - - MITK_INFO << "addPlanarFigure: size before: " << this->getNumberOfChildren(); m_PFVector->InsertElement(m_PFVector->Size(), pf); - MITK_INFO << "addPlanarFigure: size after: " << this->getNumberOfChildren(); - - - - } void mitk::PlanarFigureComposite::replaceDataNodeAt(int idx, mitk::DataNode::Pointer dn) { - MITK_INFO << "replace: size before: " << this->getNumberOfChildren(); m_DNVector->SetElement( idx, dn ); - MITK_INFO << "replace: size after: " << this->getNumberOfChildren(); - } void mitk::PlanarFigureComposite::setOperationType(PFCompositionOperation pfcOp) { this->m_compOperation = pfcOp; - MITK_INFO << "Composition set to: " << this->getOperationType(); - } mitk::PFCompositionOperation mitk::PlanarFigureComposite::getOperationType() { return this->m_compOperation; - - + + } void mitk::PlanarFigureComposite::setDisplayName(std::string displName) { m_name = displName; } std::string mitk::PlanarFigureComposite::getDisplayName() { return m_name; } int mitk::PlanarFigureComposite::getNumberOfChildren() { return m_PFVector->Size(); - + } mitk::PlanarFigure::Pointer mitk::PlanarFigureComposite::getChildAt(int idx) { - + return m_PFVector->ElementAt(idx); } mitk::DataNode::Pointer mitk::PlanarFigureComposite::getDataNodeAt(int idx) { - + return m_DNVector->ElementAt(idx); - + } //musthave implementations from superclass.... not sure if return true makes sense bool mitk::PlanarFigureComposite::SetControlPoint( unsigned int index, const Point2D &point, bool createIfDoesNotExist ) { return true; } void mitk::PlanarFigureComposite::GeneratePolyLine() { - + } void mitk::PlanarFigureComposite::GenerateHelperPolyLine(double /*mmPerDisplayUnit*/, unsigned int /*displayHeight*/) { // A circle does not require a helper object } void mitk::PlanarFigureComposite::EvaluateFeaturesInternal() { - + } void mitk::PlanarFigureComposite::PrintSelf( std::ostream& os, itk::Indent indent) const { - + } diff --git a/Modules/DiffusionImaging/IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp b/Modules/DiffusionImaging/IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp index 336a5e7066..a3f5a9bc17 100644 --- a/Modules/DiffusionImaging/IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp +++ b/Modules/DiffusionImaging/IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp @@ -1,158 +1,156 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-12-10 18:05:13 +0100 (Mi, 10 Dez 2008) $ Version: $Revision: 15922 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkNrrdQBallImageWriter.h" #include "itkMetaDataDictionary.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkImageFileWriter.h" #include "mitkImageCast.h" mitk::NrrdQBallImageWriter::NrrdQBallImageWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } mitk::NrrdQBallImageWriter::~NrrdQBallImageWriter() {} void mitk::NrrdQBallImageWriter::GenerateData() { m_Success = false; InputType* input = this->GetInput(); if (input == NULL) { itkWarningMacro(<<"Sorry, input to NrrdQBallImageWriter is NULL!"); return; } if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); io->SetFileType( itk::ImageIOBase::Binary ); io->UseCompressionOn(); typedef itk::VectorImage VecImgType; typedef itk::Image,3> ImageType; typedef itk::ImageFileWriter WriterType; WriterType::Pointer nrrdWriter = WriterType::New(); ImageType::Pointer outimage = ImageType::New(); CastToItkImage(input, outimage); VecImgType::Pointer vecImg = VecImgType::New(); vecImg->SetSpacing( outimage->GetSpacing() ); // Set the image spacing vecImg->SetOrigin( outimage->GetOrigin() ); // Set the image origin vecImg->SetDirection( outimage->GetDirection() ); // Set the image direction vecImg->SetLargestPossibleRegion( outimage->GetLargestPossibleRegion()); vecImg->SetBufferedRegion( outimage->GetLargestPossibleRegion() ); vecImg->SetVectorLength(QBALL_ODFSIZE); vecImg->Allocate(); itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); ot = ot.Begin(); itk::ImageRegionIterator it (outimage, outimage->GetLargestPossibleRegion() ); typedef ImageType::PixelType VecPixType; typedef VecImgType::PixelType VarVecType; for (it = it.Begin(); !it.IsAtEnd(); ++it) { VecPixType vec = it.Get(); VarVecType varVec(vec.GetVnlVector().data_block(), QBALL_ODFSIZE); ot.Set(varVec); ++ot; } nrrdWriter->SetInput( vecImg ); nrrdWriter->SetImageIO(io); nrrdWriter->SetFileName(m_FileName); nrrdWriter->UseCompressionOn(); try { nrrdWriter->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; } try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } m_Success = true; } void mitk::NrrdQBallImageWriter::SetInput( InputType* diffVolumes ) { this->ProcessObject::SetNthInput( 0, diffVolumes ); } mitk::QBallImage* mitk::NrrdQBallImageWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } std::vector mitk::NrrdQBallImageWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".qbi"); possibleFileExtensions.push_back(".hqbi"); return possibleFileExtensions; } diff --git a/Modules/DiffusionImaging/IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp b/Modules/DiffusionImaging/IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp index d34c76f542..ec99b1481c 100644 --- a/Modules/DiffusionImaging/IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp +++ b/Modules/DiffusionImaging/IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp @@ -1,388 +1,382 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ Version: $Revision: 18127 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef __mitkNrrdDiffusionImageReader_cpp #define __mitkNrrdDiffusionImageReader_cpp #include "mitkNrrdTbssImageReader.h" #include "itkImageFileReader.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkNiftiImageIO.h" #include #include #include #include "itksys/SystemTools.hxx" namespace mitk { template void NrrdTbssImageReader ::GenerateData() { try { // Change locale if needed const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } // READ IMAGE INFORMATION const unsigned int MINDIM = 3; const unsigned int MAXDIM = 4; MITK_INFO << "loading " << m_FileName << " via mitk::NrrdTbssImageReader... " << std::endl; // Check to see if we can read the file given the name or prefix if ( m_FileName == "" ) { itkWarningMacro( << "Filename is empty!" ) return ; } itk::NrrdImageIO::Pointer imageIO = itk::NrrdImageIO::New(); imageIO->SetFileName( m_FileName.c_str() ); imageIO->ReadImageInformation(); unsigned int ndim = imageIO->GetNumberOfDimensions(); if ( ndim < MINDIM || ndim > MAXDIM ) { itkWarningMacro( << "Sorry, only dimensions 3 is supported. The given file has " << ndim << " dimensions!" ) return; } itk::ImageIORegion ioRegion( ndim ); itk::ImageIORegion::SizeType ioSize = ioRegion.GetSize(); itk::ImageIORegion::IndexType ioStart = ioRegion.GetIndex(); unsigned int dimensions[ MAXDIM ]; dimensions[ 0 ] = 0; dimensions[ 1 ] = 0; dimensions[ 2 ] = 0; dimensions[ 3 ] = 0; float spacing[ MAXDIM ]; spacing[ 0 ] = 1.0f; spacing[ 1 ] = 1.0f; spacing[ 2 ] = 1.0f; spacing[ 3 ] = 1.0f; Point3D origin; origin.Fill(0); unsigned int i; for ( i = 0; i < ndim ; ++i ) { ioStart[ i ] = 0; ioSize[ i ] = imageIO->GetDimensions( i ); if(iGetDimensions( i ); spacing[ i ] = imageIO->GetSpacing( i ); if(spacing[ i ] <= 0) spacing[ i ] = 1.0f; } if(i<3) { origin[ i ] = imageIO->GetOrigin( i ); } } ioRegion.SetSize( ioSize ); ioRegion.SetIndex( ioStart ); - MITK_INFO << "ioRegion: " << ioRegion << std::endl; imageIO->SetIORegion( ioRegion ); void* buffer = new unsigned char[imageIO->GetImageSizeInBytes()]; imageIO->Read( buffer ); //mitk::Image::Pointer static_cast(this->GetOutput())image = mitk::Image::New(); if((ndim==4) && (dimensions[3]<=1)) ndim = 3; if((ndim==3) && (dimensions[2]<=1)) ndim = 2; mitk::PixelType pixelType( imageIO->GetComponentTypeInfo(), imageIO->GetNumberOfComponents(), imageIO->GetPixelType() ); static_cast(this->GetOutput())->Initialize( pixelType, ndim, dimensions ); static_cast(this->GetOutput())->SetImportChannel( buffer, 0, Image::ManageMemory ); // access direction of itk::Image and include spacing mitk::Matrix3D matrix; matrix.SetIdentity(); unsigned int j, itkDimMax3 = (ndim >= 3? 3 : ndim); for ( i=0; i < itkDimMax3; ++i) for( j=0; j < itkDimMax3; ++j ) matrix[i][j] = imageIO->GetDirection(j)[i]; // re-initialize PlaneGeometry with origin and direction PlaneGeometry* planeGeometry = static_cast (static_cast (this->GetOutput())->GetSlicedGeometry(0)->GetGeometry2D(0)); planeGeometry->SetOrigin(origin); planeGeometry->GetIndexToWorldTransform()->SetMatrix(matrix); // re-initialize SlicedGeometry3D SlicedGeometry3D* slicedGeometry = static_cast(this->GetOutput())->GetSlicedGeometry(0); slicedGeometry->InitializeEvenlySpaced(planeGeometry, static_cast(this->GetOutput())->GetDimension(2)); slicedGeometry->SetSpacing(spacing); // re-initialize TimeSlicedGeometry static_cast(this->GetOutput())->GetTimeSlicedGeometry()->InitializeEvenlyTimed(slicedGeometry, static_cast(this->GetOutput())->GetDimension(3)); buffer = NULL; - MITK_INFO << "number of image components: "<< static_cast(this->GetOutput())->GetPixelType().GetNumberOfComponents() << std::endl; - - + //MITK_INFO << "number of image components: "<< static_cast(this->GetOutput())->GetPixelType().GetNumberOfComponents() << std::endl; // READ TBSS HEADER INFORMATION typename ImageType::Pointer img; std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".tbss") { typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(this->m_FileName); reader->SetImageIO(imageIO); reader->Update(); 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("tbss") != std::string::npos) { - MITK_INFO << *itKey << " ---> " << metaString; + //MITK_INFO << *itKey << " ---> " << metaString; if(metaString == "ROI") { - MITK_INFO << "Read the ROI info"; + //MITK_INFO << "Read the ROI info"; ReadRoiInfo(imgMetaDictionary); // move back into if statement static_cast(this->GetOutput())->SetTbssType(OutputType::ROI); } } } } // RESET LOCALE try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } - MITK_INFO << "...finished!" << std::endl; + //MITK_INFO << "...finished!" << std::endl; } catch(std::exception& e) { - MITK_INFO << "Std::Exception while reading file!!"; - MITK_INFO << e.what(); + //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 void NrrdTbssImageReader ::ReadRoiInfo(itk::MetaDataDictionary dict) { std::vector imgMetaKeys = dict.GetKeys(); std::vector::const_iterator itKey = imgMetaKeys.begin(); std::string metaString; std::vector< itk::Index<3> > roi; for (; itKey != imgMetaKeys.end(); itKey ++) { double x,y,z; itk::Index<3> ix; itk::ExposeMetaData (dict, *itKey, metaString); if (itKey->find("ROI_index") != std::string::npos) { - MITK_INFO << *itKey << " ---> " << metaString; + //MITK_INFO << *itKey << " ---> " << metaString; sscanf(metaString.c_str(), "%lf %lf %lf\n", &x, &y, &z); ix[0] = x; ix[1] = y; ix[2] = z; roi.push_back(ix); } else if(itKey->find("preprocessed FA") != std::string::npos) { - MITK_INFO << *itKey << " ---> " << metaString; + //MITK_INFO << *itKey << " ---> " << metaString; static_cast(this->GetOutput())->SetPreprocessedFA(true); static_cast(this->GetOutput())->SetPreprocessedFAFile(metaString); } // Name of structure if (itKey->find("structure") != std::string::npos) { - MITK_INFO << *itKey << " ---> " << metaString; + //MITK_INFO << *itKey << " ---> " << metaString; static_cast(this->GetOutput())->SetStructure(metaString); } } static_cast(this->GetOutput())->SetRoi(roi); } template const char* NrrdTbssImageReader ::GetFileName() const { return m_FileName.c_str(); } template void NrrdTbssImageReader ::SetFileName(const char* aFileName) { m_FileName = aFileName; } template const char* NrrdTbssImageReader ::GetFilePrefix() const { return m_FilePrefix.c_str(); } template void NrrdTbssImageReader ::SetFilePrefix(const char* aFilePrefix) { m_FilePrefix = aFilePrefix; } template const char* NrrdTbssImageReader ::GetFilePattern() const { return m_FilePattern.c_str(); } template void NrrdTbssImageReader ::SetFilePattern(const char* aFilePattern) { m_FilePattern = aFilePattern; } template bool NrrdTbssImageReader ::CanReadFile(const std::string filename, const std::string filePrefix, const std::string filePattern) { // First check the extension if( filename == "" ) return false; // check if image is serie if( filePattern != "" && filePrefix != "" ) return false; std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".tbss") { 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(); return false; } typename ImageType::Pointer img = reader->GetOutput(); - itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); + 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("tbss") != std::string::npos) { if (metaString.find("ROI") != std::string::npos) { return true; } } } } return false; } } //namespace MITK #endif diff --git a/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp b/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp index b5f57ff7cb..e30b2a504b 100644 --- a/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp +++ b/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp @@ -1,280 +1,275 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ Version: $Revision: 18127 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkNrrdTensorImageReader.h" #include "itkImageFileReader.h" #include "itkImageRegionIterator.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkDiffusionTensor3D.h" #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 { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } 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 ++) { itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); if (itKey->find("measurement frame") != std::string::npos) { - MITK_INFO << *itKey << " ---> " << metaString.c_str(); 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(); - - MITK_INFO << "Will apply following measurement frame: \n" << measFrame; } } } if (numComponents==6) { while (!it.IsAtEnd()) { // 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 = tensor.PreMultiply(measFrame); tensor = tensor.PostMultiply(measFrameTransp); fixVec = tensor; } ot.Set(fixVec); ++ot; ++it; } } else if(numComponents==9) { while (!it.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) { tensor = tensor.PreMultiply(measFrame); tensor = tensor.PostMultiply(measFrameTransp); } FixPixType fixVec(tensor); ot.Set(fixVec); ++ot; ++it; } } else { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!"); } this->GetOutput()->InitializeByItk(vecImg.GetPointer()); this->GetOutput()->SetVolume(vecImg->GetBufferPointer()); try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; 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; } } //namespace MITK diff --git a/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp b/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp index 4a9bc8dd10..3ab6edd85c 100644 --- a/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp +++ b/Modules/DiffusionImaging/IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp @@ -1,130 +1,128 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-12-10 18:05:13 +0100 (Mi, 10 Dez 2008) $ Version: $Revision: 15922 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkNrrdTensorImageWriter.h" #include "itkMetaDataDictionary.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkImageFileWriter.h" #include "itkDiffusionTensor3D.h" #include "mitkImageCast.h" mitk::NrrdTensorImageWriter::NrrdTensorImageWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } mitk::NrrdTensorImageWriter::~NrrdTensorImageWriter() {} void mitk::NrrdTensorImageWriter::GenerateData() { m_Success = false; InputType* input = this->GetInput(); if (input == NULL) { itkWarningMacro(<<"Sorry, input to NrrdTensorImageWriter is NULL!"); return; } if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); if ( locale.compare(currLocale)!=0 ) { try { - MITK_INFO << " ** Changing locale from " << setlocale(LC_ALL, NULL) << " to '" << locale << "'"; setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); io->SetFileType( itk::ImageIOBase::Binary ); io->UseCompressionOn(); typedef itk::Image,3> ImageType; typedef itk::ImageFileWriter WriterType; WriterType::Pointer nrrdWriter = WriterType::New(); ImageType::Pointer outimage = ImageType::New(); CastToItkImage(input, outimage); nrrdWriter->SetInput( outimage ); nrrdWriter->SetImageIO(io); nrrdWriter->SetFileName(m_FileName); nrrdWriter->UseCompressionOn(); try { nrrdWriter->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; } try { - MITK_INFO << " ** Changing locale back from " << setlocale(LC_ALL, NULL) << " to '" << currLocale << "'"; setlocale(LC_ALL, currLocale.c_str()); } catch(...) { MITK_INFO << "Could not reset locale " << currLocale; } m_Success = true; } void mitk::NrrdTensorImageWriter::SetInput( InputType* diffVolumes ) { this->ProcessObject::SetNthInput( 0, diffVolumes ); } mitk::TensorImage* mitk::NrrdTensorImageWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } std::vector mitk::NrrdTensorImageWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".dti"); possibleFileExtensions.push_back(".hdti"); return possibleFileExtensions; } diff --git a/Modules/DiffusionImaging/Rendering/mitkDiffusionImageMapper.cpp b/Modules/DiffusionImaging/Rendering/mitkDiffusionImageMapper.cpp index bd95ac041f..365b011675 100644 --- a/Modules/DiffusionImaging/Rendering/mitkDiffusionImageMapper.cpp +++ b/Modules/DiffusionImaging/Rendering/mitkDiffusionImageMapper.cpp @@ -1,67 +1,63 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-05-12 19:56:03 +0200 (Di, 12 Mai 2009) $ Version: $Revision: 17179 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef DiffusionImageMapper_txx_HEADER_INCLUDED #define DiffusionImageMapper_txx_HEADER_INCLUDED #include "mitkProperties.h" #include "mitkDiffusionImage.h" template mitk::DiffusionImageMapper::DiffusionImageMapper() { - MITK_INFO << "DiffMapper init"; + } template mitk::DiffusionImageMapper::~DiffusionImageMapper() { - MITK_INFO << "DiffMapper destroyed"; + } template void mitk::DiffusionImageMapper::GenerateDataForRenderer( mitk::BaseRenderer *renderer ) { int displayIndex(0); this->GetDataNode()->GetIntProperty( "DisplayChannel", displayIndex, renderer ); mitk::Image *input = const_cast< mitk::Image* >( this->GetInput() ); mitk::DiffusionImage *input2 = dynamic_cast< mitk::DiffusionImage* >( input ); - MITK_INFO << "displayindex: " << displayIndex; - - - input2->SetDisplayIndexForRendering(displayIndex); Superclass::GenerateDataForRenderer(renderer); } template void mitk::DiffusionImageMapper::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite) { node->AddProperty( "DisplayChannel", mitk::IntProperty::New( 0 ), renderer, overwrite ); Superclass::SetDefaultProperties(node, renderer, overwrite); } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/reparametrize_arclen2.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/reparametrize_arclen2.cpp index 657f72491b..d9f5c2cf06 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/reparametrize_arclen2.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/reparametrize_arclen2.cpp @@ -1,97 +1,91 @@ #include "mitkParticle.h" #include "auxilary_classes.cpp" using namespace std; #define PI 3.1415926536 float squareNorm( vnl_vector_fixed v ) { return v[0]*v[0]+v[1]*v[1]+v[2]*v[2]; } mitk::Particle* createMitkParticle(Particle* p) { mitk::Particle* particle = new mitk::Particle(); pVector R = p->R; vnl_vector_fixed pos; pos[0] = R[0]-0.5; pos[1] = R[1]-0.5; pos[2] = R[2]-0.5; particle->SetPosition(pos); - //MITK_INFO << "mitk: " << particle->GetPosition(); pVector N = p->N; vnl_vector_fixed dir; dir[0] = N[0]; dir[1] = N[1]; dir[2] = N[2]; particle->SetDirection(dir); return particle; } Particle* createParticle(Particle* particle){ Particle* p = new Particle(); p->R = pVector(particle->R); p->N = pVector(particle->N); return p; } vector* ResampleFibers(vector* particleContainer, float len) { //vector< mitk::Particle::Pointer >* outContainer = new::vector< mitk::Particle::Pointer >(); vector< Particle* >* container = new vector< Particle* >(); int numPoints = particleContainer->size(); if(numPoints<2) return container; Particle* source = createParticle(particleContainer->at(0)); Particle* sink = createParticle(particleContainer->at(numPoints-1)); float Leng = 0; container->push_back(source); - // source->R.storeXYZ(); - // MITK_INFO << source->R.store[0] << ", " << source->R.store[1] << ", " << source->R.store[2]; float dtau = 0; int cur_p = 1; int cur_i = 1; pVector dR; float normdR; for (;;) { while (dtau <= len && cur_p < numPoints) { dR = particleContainer->at(cur_p)->R - particleContainer->at(cur_p-1)->R; normdR = sqrt(dR.norm_square()); dtau += normdR; Leng += normdR; cur_p++; } if (dtau >= len) // if particles reach next voxel { - // proposal = current particle position - (current p - last p)*(current fibre length - len???)/(norm current p-last p) Particle* p = new Particle(); p->R = particleContainer->at(cur_p-1)->R - dR*( (dtau-len)/normdR ); p->N = dR; - // p->R.storeXYZ(); - // MITK_INFO << p->R.store[0] << ", " << p->R.store[1] << ", " << p->R.store[2]; container->push_back(p); } else { container->push_back(sink); break; } dtau = dtau-len; cur_i++; } return container; } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index 290a6243fd..c95dec827f 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,485 +1,485 @@ #include "itkGibbsTrackingFilter.h" #include #include #include "itkPointShell.h" #include "GibbsTracking/BuildFibres.cpp" #pragma GCC visibility push(default) #include #pragma GCC visibility pop #include #include #include #include #include #include #include "GibbsTracking/reparametrize_arclen2.cpp" #include #include struct LessDereference { template bool operator()(const T * lhs, const T * rhs) const { return *lhs < *rhs; } }; namespace itk{ template< class TInputOdfImage, class TInputROIImage > GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GibbsTrackingFilter(): m_TempStart(0.1), m_TempEnd(0.001), m_NumIt(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ChempotConnection(10), m_ChempotParticle(0), m_InexBalance(0), m_Chempot2(0.2), m_FiberLength(10), m_AbortTracking(false), m_NumConnections(0), m_NumParticles(0), m_NumAcceptedFibers(0), m_CurrentStep(0), m_SubtractMean(true), m_BuildFibers(false), m_Sampler(NULL), m_Steps(10), m_Memory(0), m_ProposalAcceptance(0) { //this->m_MeasurementFrame.set_identity(); this->SetNumberOfRequiredInputs(2); //Filter needs a DWI image + a Mask Image } template< class TInputOdfImage, class TInputROIImage > GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::~GibbsTrackingFilter(){ delete BESSEL_APPROXCOEFF; if (m_Sampler!=NULL) delete m_Sampler; } template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::ComputeFiberCorrelation(){ // float bD = 15; // vnl_matrix_fixed bDir = // *itk::PointShell >::DistributePointShell(); // const int N = QBALL_ODFSIZE; // vnl_matrix_fixed temp = bDir.transpose(); // vnl_matrix_fixed C = temp*bDir; // vnl_matrix_fixed Q = C; // vnl_vector_fixed mean; // for(int i=0; i repMean; // for (int i=0; i P = Q*Q; // std::vector pointer; // pointer.reserve(N*N); // double * start = C.data_block(); // double * end = start + N*N; // for (double * iter = start; iter != end; ++iter) // { // pointer.push_back(iter); // } // std::sort(pointer.begin(), pointer.end(), LessDereference()); // vnl_vector_fixed alpha; // vnl_vector_fixed beta; // for (int i=0; im_Meanval_sq = (sum*sum)/N; // vnl_vector_fixed alpha_0; // vnl_vector_fixed alpha_2; // vnl_vector_fixed alpha_4; // vnl_vector_fixed alpha_6; // for(int i=0; i T; // T.set_column(0,alpha_0); // T.set_column(1,alpha_2); // T.set_column(2,alpha_4); // T.set_column(3,alpha_6); // vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; -// MITK_INFO << "Bessel oefficients: " << coeff; +// MITK_INFO << "itkGibbsTrackingFilter: Bessel oefficients: " << coeff; BESSEL_APPROXCOEFF = new float[4]; // BESSEL_APPROXCOEFF[0] = coeff(0); // BESSEL_APPROXCOEFF[1] = coeff(1); // BESSEL_APPROXCOEFF[2] = coeff(2); // BESSEL_APPROXCOEFF[3] = coeff(3); BESSEL_APPROXCOEFF[0] = -0.1714; BESSEL_APPROXCOEFF[1] = 0.5332; BESSEL_APPROXCOEFF[2] = -1.4889; BESSEL_APPROXCOEFF[3] = 2.0389; } // build fibers from tracking result template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::BuildFibers(float* points, int numPoints) { - MITK_INFO << "Building fibers ..."; + MITK_INFO << "itkGibbsTrackingFilter: Building fibers ..."; typename InputQBallImageType::Pointer odfImage = dynamic_cast(this->GetInput(0)); double spacing[3]; spacing[0] = odfImage->GetSpacing().GetElement(0); spacing[1] = odfImage->GetSpacing().GetElement(1); spacing[2] = odfImage->GetSpacing().GetElement(2); // initialize array of particles CCAnalysis ccana(points, numPoints, spacing); // label the particles according to fiber affiliation and return number of fibers int numFibers = ccana.iterate(m_FiberLength); if (numFibers<=0){ - MITK_INFO << "0 fibers accepted"; + MITK_INFO << "itkGibbsTrackingFilter: 0 fibers accepted"; return; } // fill output datastructure m_FiberBundle.clear(); for (int i = 0; i < numFibers; i++) { vector< Particle* >* particleContainer = ccana.m_FiberContainer->at(i); // resample fibers std::vector< Particle* >* pCon = ResampleFibers(particleContainer, 0.9*spacing[0]); FiberTractType tract; for (int j=0; jsize(); j++) { Particle* particle = pCon->at(j); pVector p = particle->R; itk::Point point; point[0] = p[0]-0.5; point[1] = p[1]-0.5; point[2] = p[2]-0.5; tract.push_back(point); delete(particle); } m_FiberBundle.push_back(tract); delete(pCon); } m_NumAcceptedFibers = numFibers; MITK_INFO << "itkGibbsTrackingFilter: " << numFibers << " fibers accepted"; } // fill output fiber bundle datastructure template< class TInputOdfImage, class TInputROIImage > typename GibbsTrackingFilter< TInputOdfImage, TInputROIImage >::FiberBundleType* GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return &m_FiberBundle; } // get memory allocated for particle grid template< class TInputOdfImage, class TInputROIImage > float GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GetMemoryUsage() { if (m_Sampler!=NULL) return m_Sampler->m_ParticleGrid.GetMemoryUsage(); return 0; } // perform global tracking template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GenerateData(){ // input qball image m_ItkQBallImage = dynamic_cast(this->GetInput(0)); // approximationscoeffizienten der // teilchenkorrelationen im orientierungsraum // 4er vektor ComputeFiberCorrelation(); // image sizes and spacing int qBallImageSize[4] = {QBALL_ODFSIZE, m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(0), m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(1), m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(2)}; double qBallImageSpacing[3] = {m_ItkQBallImage->GetSpacing().GetElement(0),m_ItkQBallImage->GetSpacing().GetElement(1),m_ItkQBallImage->GetSpacing().GetElement(2)}; // make sure image has enough slices if (qBallImageSize[1]<3 || qBallImageSize[2]<3 || qBallImageSize[3]<3) { - MITK_INFO << "image size < 3 not supported"; + MITK_INFO << "itkGibbsTrackingFilter: image size < 3 not supported"; m_AbortTracking = true; } // calculate rotation matrix vnl_matrix_fixed directionMatrix = m_ItkQBallImage->GetDirection().GetVnlMatrix(); vnl_vector_fixed d0 = directionMatrix.get_column(0); d0.normalize(); vnl_vector_fixed d1 = directionMatrix.get_column(1); d1.normalize(); vnl_vector_fixed d2 = directionMatrix.get_column(2); d2.normalize(); directionMatrix.set_column(0, d0); directionMatrix.set_column(1, d1); directionMatrix.set_column(2, d2); vnl_matrix_fixed I = directionMatrix*directionMatrix.transpose(); if(!I.is_identity(mitk::eps)){ - MITK_INFO << "Image direction is not a rotation matrix. Tracking not possible!"; + MITK_INFO << "itkGibbsTrackingFilter: image direction is not a rotation matrix. Tracking not possible!"; m_AbortTracking = true; } // generate local working copy of image buffer int bufferSize = qBallImageSize[0]*qBallImageSize[1]*qBallImageSize[2]*qBallImageSize[3]; float* qBallImageBuffer = (float*) m_ItkQBallImage->GetBufferPointer(); float* workingQballImage = new float[bufferSize]; for (int i=0; i0 && i%qBallImageSize[0] == 0 && i>0) { sum /= qBallImageSize[0]; for (int j=i-qBallImageSize[0]; jGetBufferPointer(); maskImageSize[0] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(0); maskImageSize[1] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(1); maskImageSize[2] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(2); } else { mask = 0; maskImageSize[0] = qBallImageSize[1]; maskImageSize[1] = qBallImageSize[2]; maskImageSize[2] = qBallImageSize[3]; } int mask_oversamp_mult = maskImageSize[0]/qBallImageSize[1]; // load lookuptable QString applicationDir = QCoreApplication::applicationDirPath(); if (applicationDir.endsWith("bin")) applicationDir.append("/"); else applicationDir.append("\\..\\"); ifstream BaryCoords; QString lutPath(applicationDir+"FiberTrackingLUTBaryCoords.bin"); BaryCoords.open(lutPath.toStdString().c_str(), ios::in | ios::binary); float* coords; if (BaryCoords.is_open()) { float tmp; coords = new float [1630818]; BaryCoords.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { BaryCoords.read((char *)&tmp, sizeof(tmp)); coords[i] = tmp; } BaryCoords.close(); } else { - MITK_INFO << "Unable to open barycoords file"; + MITK_INFO << "itkGibbsTrackingFilter: unable to open barycoords file"; m_AbortTracking = true; } ifstream Indices; lutPath = applicationDir+"FiberTrackingLUTIndices.bin"; Indices.open(lutPath.toStdString().c_str(), ios::in | ios::binary); int* ind; if (Indices.is_open()) { int tmp; ind = new int [1630818]; Indices.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { Indices.read((char *)&tmp, 4); ind[i] = tmp; } Indices.close(); } else { - MITK_INFO << "Unable to open indices file"; + MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; m_AbortTracking = true; } // initialize sphere interpolator with lookuptables SphereInterpolator *sinterp = new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); // get paramters float minSpacing; if(qBallImageSpacing[0]m_NumIt) { - MITK_INFO << "not enough iterations!"; + MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); // setup metropolis hastings sampler MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; if (m_Sampler!=NULL) delete m_Sampler; m_Sampler = new RJMCMC(NULL, 0, workingQballImage, qBallImageSize, qBallImageSpacing, cellsize); // setup energy computer MITK_INFO << "itkGibbsTrackingFilter: setting up Energy-computer"; EnergyComputer encomp(workingQballImage,qBallImageSize,qBallImageSpacing,sinterp,&(m_Sampler->m_ParticleGrid),mask,mask_oversamp_mult, directionMatrix); encomp.setParameters(m_ParticleWeight,m_ParticleWidth,m_ChempotConnection*m_ParticleLength*m_ParticleLength,m_ParticleLength,curvatureHardThreshold,m_InexBalance,m_Chempot2); m_Sampler->SetEnergyComputer(&encomp); m_Sampler->SetParameters(m_TempStart,singleIts,m_ParticleLength,curvatureHardThreshold,m_ChempotParticle); // main loop for( int step = 0; step < m_Steps; step++ ) { if (m_AbortTracking) break; m_CurrentStep = step+1; float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); - MITK_INFO << "iterating step " << m_CurrentStep; + MITK_INFO << "itkGibbsTrackingFilter: iterating step " << m_CurrentStep; m_Sampler->SetTemperature(temperature); m_Sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); - MITK_INFO << "proposal acceptance: " << 100*m_ProposalAcceptance << "%"; - MITK_INFO << "particles: " << m_NumParticles; - MITK_INFO << "connections: " << m_NumConnections; - MITK_INFO << "progress: " << 100*(float)step/m_Steps << "%"; + MITK_INFO << "itkGibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; + MITK_INFO << "itkGibbsTrackingFilter: particles: " << m_NumParticles; + MITK_INFO << "itkGibbsTrackingFilter: connections: " << m_NumConnections; + MITK_INFO << "itkGibbsTrackingFilter: progress: " << 100*(float)step/m_Steps << "%"; if (m_BuildFibers) { int numPoints = m_Sampler->m_ParticleGrid.pcnt; float* points = new float[numPoints*m_Sampler->m_NumAttributes]; m_Sampler->WriteOutParticles(points); BuildFibers(points, numPoints); delete points; m_BuildFibers = false; } } int numPoints = m_Sampler->m_ParticleGrid.pcnt; float* points = new float[numPoints*m_Sampler->m_NumAttributes]; m_Sampler->WriteOutParticles(points); BuildFibers(points, numPoints); delete points; delete sinterp; delete coords; delete ind; delete workingQballImage; m_AbortTracking = true; m_BuildFibers = false; - MITK_INFO << "done generate data"; + MITK_INFO << "itkGibbsTrackingFilter: done generate data"; } }